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ABSTRACT 



We investigate the evolution of galactic discs in N-body Tree-SPH simulations. We find that discs, initially truncated at three scale- 
lengths, can triple their radial extent, solely driven by secular evolution. At the same time, the initial radial metallicity gradients are 
flattened and even reversed in the outer discs. Both Type I (single exponential) and Type 11 (down-turning) observed disc surface- 
brightness profiles can be explained by our findings. We show that profiles with breaks beyond the bar's outer Lindblad resonance, at 
present only explained as the effect of star-formation threshold, can occur even if no star formation is considered. We explain these 
results with the strong angular momentum outward transfer, resulting from torques and radial migration associated with multiple 
patterns, such as central bars and spiral waves of different multiplicity. We find that even for stars ending up on cold orbits, the changes 
in angular momentum exhibit complex structure as a function of radius, unlike the expected effect of transient spirals alone. We show 
that the bars in all of our simulations are the most effective drivers of radial migration through their corotation resonance, throughout 
the 3 Gyr of evolution studied. Focussing on one of our models, we find evidence for non-linear coupling among m = 1, 2, 3 and 4 
density waves, where m is the pattern multiplicity. In this way the waves involved conspire to carry the energy and angular momentum 
extracted by the first mode from the inner parts of the disc much farther out than a single mode could. We suggest that the naturally 
occurring larger resonance widths at galactic radii beyond four scale-lengths may have profound consequences on the formation and 
location of breaks in disc density profiles, provided spirals are present at such large distances. We also consider the effect of gas 
inflow and show that when in-plane smooth gas accretion of ~ 5 Mo/yr is included, the outer discs become more unstable, leading to a 
strong increase in the stellar velocity dispersion. This, in turn, causes the formation of a Type 111 (up-turning) profile in the old stellar 
population. We propose that observations of Type 111 surface brightness profiles, combined with an up-tum in the stellar velocity 
dispersions beyond the disc break, could be a signature of ongoing gas-accretion. The results of this study suggest that disc outskirts 
comprised of stars migrated from the inner disc would have relatively large radial velocity dispersions (> 30 km/s at 6 scale-lengths 
for Milky Way-size systems), and significant thickness when seen edge-on. 

Key words. Keywords should be given 



1. Introduction 

It has been recognized for a long time now that the ra- 
dial distribution of the light of stellar discs of most spiral 
galaxies i s wel l approximated by an ex ponential profi l e (e.g . 
iFree man; 19 70|). The classical wo rks of Ivan der KruitI (1 19791) 
and Ivan der Kruit and Searld (|1981|) have shown that some spi- 
ral, mostly late-type edge-on galaxies, have truncated discs (i.e., 
the surface-brightness profile shows a very sharp exponential 
fall-off beyond the truncation radius, typically 2-4 exponential 
scale-lengths). A number of recent studies looking into galax- 
ies at a range of inclinations have actually found that galac- 
tic discs do not end at the truncation, but rather exhibit a 
change in the exponential scale-length of the light distribution 
(e.g. Pohlen et al., 2002; Erwin et al. , 2005 ; Pohlen and Tmiill o, 
120061; iHunter and Elmegreeni l2006b . It is now acknowledged 
that the surface-brightness profiles of the mass-carrying old 



stellar population in the outer discs can be a simple continu- 
ation of the inner exponential profiles (Type I), steeper (Type 
II) or even shallower (Type III) than the extrapolations of the 
inner profiles to large radii (se e JBland-Hawthorn et al,i l2005t 
lEllis and Bland-Hawthornl2007l and lErwin et al.ll2008l hereafter 
EPB08, and references therein). 

The diversity of disc surface-brightness profiles is observed 
for both early-type and late-type disc galaxies. The frequency of 
different shapes of surface-brightness profiles appears to corre- 
late, however, with galaxy morphological type. Whereas discs 
for which the outer surface-brightness profiles are steeper than 
the extrapolation of the inner profiles (Type II) are more frequent 
in late-type galaxies, the fraction of galaxies with shallower 
outer profiles (Type III) rises toward earlier morphological types 
(Pohle n and Tmiillo 2006, EPB08). Galactic discs with single 
exponential profiles (Type I) appear to be much rarer in late-type 
galaxies (EPB08). For the classical truncated surface-brightness 
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profiles the breaks occur at small er normalized radii in late -t ype 
galaxies than early-type galaxies (iPohlen and Truiilloll2006t) . 

By extending t he EPB08 samp l e of 66 barred galaxies with 
47 unbarred ones, [Gutierrez et al.l (1201 ll ) made the first clear 
statements about the trends of outer-disc profile types along the 
Hubble sequence and their global frequencies. They found that 
Type I profiles are most frequently found in early-type discs, de- 
creasing from one-third of all SO-Sa discs to 10% of the latest 
type spirals. Conversely, Type II profiles increase in frequency 
with Hubble type, from ~ 25% of SO galaxies to ~ 80% of Sd- 
Sm spirals. The fractions of Type I, II, and III profiles for all disc 
galaxies (Hubble types SO-Sm) were found to be 21%, 50%, and 
38%, respectively. 

Cosmic environment could also play a n important role in 
shaping disc profiles. In a recent study, lErwin et all (1201 2|) 
showed that while SO field galaxies present statistically almost 
equal shares of Types I, II, and III, the Type II discs are missing 
for a sample selected from the Virgo cluster If this were caused 
by environmental effects, the authors proposed two solutions: ei- 
ther something in the cluster (or protocluster) environment trans- 
forms Type II profiles into Type I, or else something prevents a 
Type I to Type II transition which is common in the field. 

Resolving the stellar contents of galaxies gives access 
to much fainter surface-brightness levels and allows for 
the estimati on of stellar po pulation ages and metallicities 
dBland-Hawthorn et al.Ll20()5l) . This could provide an insight on 
the processes giving rise to the properties of the outskirts of 
galactic discs . In the ir study of the low mass galaxy NGC 4244, 
Ide Jonp et alj (l2007h have shown that the break in the stellar 
count profiles occurs at the same radius for young, interme- 
diate age, and old stars, independently from the height above 
the galactic plane. The low-mass disc galaxy M33 has a trun- 
cated disc (iFerguson et al.L l2007h . Measurements of the radial 
change of stellar ages of the dominant stellar populations show 
that the inner negative stellar age gradient reverses at the ra- 
dius where the break of the surface-brightness p rofile is mea- 
sured (iBarker et al.Ll2007bHwrmams et al.ll2009l) . The metalUc- 
ity gradient beyond the truncation radius is consistent however 
with the extrapolation of the metallicity gradient of the inner 
disc (iBarker et al.L l2007ah . In NGC 300, a galaxy with compa- 
rable properties to those of M33, a single exponential surface- 
brightness profile extends out to ten scale lengths of the galaxy 
(iBland-Hawthorn et al1.l2005h . Similarly to M33, the inner neg- 
ative metallicity gradient of NGC 300 is found to extend out to 
6-7 s cale lengths, but b eyond whic h the metallicity gr adient flat- 
tens (IVlaiic et al.l2009l and see also lVlaiic et al.l20I ll for similar 
results in NGC 7793). 

A number of scenarios have been investigated to explain the 
observed diversity of surface-brightness profiles seen in galactic 
discs, which can be separated into two categories. The first cat- 
egory advocates that the light profiles result from stars formed 
in the observed distribution, as a resu lts of either a limit in the 
gas distribution (Ivan de r Kruit. 1987), o r by imposing a thresh- 
old for star formation (IKennicutt .19891) . lElmegreen and Hunted 
(|2006|) have shown that a double-exponential profile could result 
from a multicomponent star formation prescription, where tur- 
bulent compression in the outskiits of the disc allows for cloud 
formation and star formation despite subcritical densities. 

The second category of models speculate that stars in the 
outer d iscs were deposited in those regions by dynamical pro- 
cesses. ISellwood and BinnevI (l2002h have invoked the idea of 
radial migration from transient spiral arms, where a large num- 
ber of transients can mix stars at the corotation radius (CR) ra- 
dially with the largest migration coinciding with the growth of 



strong spirals. Stars beyond the break in the exponential pro- 
file can thus be scattered outward from the inner disc. Since 
old stars have longer timescales to populate the outer regions 
of galactic discs, this leads to a change in the radial mean stel- 
lar age profile at the break radius and a positive age gradient 
in outer disc (iRoskar et al.L |2008|) . The latter theoretical pre- 
dic tions have found good agreement with observations (e.g., 
Yoachim et al.l l2010i |2012|) . By analy zing high-resolution hy- 
drodyn amical disc galaxy simulations ISanchez-Blazquez et al.l 
(l2009h have argued, however, that the change in the age gra- 
dient at the break radius may, alternatively, be due to a different 
star formation rate between regions inside the break radius and 
those outside. The efficiency of radial migration caused by tran- 
sient spiral density waves is suspected to be significantly lower 
in low-mass discs than in Milky Way (MW)-like galaxies, as 
stron g transient spirals are not expected to grow in such sys- 
tems dGogarten et al.L |2010|) . More recent comparison between 
the migra tion seen in MW-size and lo w-mass systems was pre- 
sented bv lRadburn-Smith et al.l (l2012h . where observational evi- 
dence of radial migration was found for NGC 7793. 

Minor mergers are also expected to play a role in populat- 
the outer regions of galactic discs (IPenarrubia et al I l2006t 
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lOuiUen et al.Ll2009HBird et al.Ll2012HPurceII et aI.Ll20Tll) . In the 
MW, small satelli tes on radial, i n-plan e orbits can cause mixing 
in the outer disc (lOuillen et al.L 12009') and thus account for the 
fracti on of low-metalli city stars present in the solar neighbor- 
hood (iHavwoodL l2008h . Such perturbations also create transient 
spiral s, giving rise to structure in the velocity (iMinchev et al. 



2009h and the en ergy-angular momentu m space (iGomez et al. 



Recently, iGomez et al.l (1201 2bl) related features in the 
energy distribution of the SEGUE G-dwarf stellar sample to 
disc perturbations by the Sagittarius dwarf galaxy, also show- 
ing that this MW satellite may be responsible for vertical den- 
sity waves near the Sun Gomez et al. (2012c). Given that the ef- 
fect of minor mergers can be seen near the Sun, perturbations on 
the outer discs mu s t be e ven larger due to the lower disc den- 
sity. lYounger et al.l (l2007l) have demonstrated that the observed 
excess surface-brightness at large radii relative to the inner ex- 
ponential profile in some galaxies could be the result of minor 
mergers. This requires, however, a significant gas supply in the 
primary disc and the encounter to be prograde with moderate 
orbital angular momentum. 

The increase of the fraction of galactic discs with antitrun- 
cated profiles toward earlier morphological types (EPB08), and 
the t endency of bar s to be larger in early-typ e disc galaxies 
(e.g. IChapelon et"an I1999L lErwin et al.L 120051) could be used 
to argue that the presence of a bar might play an important 
role in shaping the propertie s of t he outskirts of g alactic discs. 
iMinchev and FamaevI (1201 Ol) and iMinchev et al.l (1201 lal) have 
shown that spiral structure interacting with a central bar can be 
an extremely efficient mechanism for radial mixing in galactic 
discs. Although angular momentum changes are most prominent 
in the vici nity of the CR radius of eac h individual perturber(as 
shown bv lSellwood and Binnevll2002l) . the non-linear coupling 
between the bar and spiral waves (e.g., iTagger et al.l 119871 
ISvgnet et al.lll988l) makes this mechanism effective over the en- 
tire galactic disc. Being non-linear, this new way of mixing can 
be significantly more efficient at increasing the angular momen- 
tum than transient sp irals alone and works with bo th short- and 
long-lived spirals (see lMinchev and Famaevir20I0L for a detailed 
discussion). Subsequent studies of radial diffusion coefficients 
have c onfirmed the validity of bar-spira l coupling dShevchenka 
I2OIIL iBrunettietaLL 1201 ll) . Recendy, IComparetta and OuillenI 
(I2OI2I) showed that, even if patterns are long-lived, radial mi- 
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Table 1. Galaxy modeling parameters 
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Mh [2.3 X 10' Mo] 
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0.1 


0.2 
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rg [kpc] 
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rn [kpc] 




10 


10 


12 


10 


12 


a., [kpc] 
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b. [kpc] 




0.5 
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0.5 


0.5 


0.5 


flg [kpc] 




- 


5 


6 


5 


6 


bg [kpc] 




- 


0.2 


0.2 


0.2 


0.2 


N, 




- 


80000 


160 000 


80000 


80000 


N, 




320 000 


240 000 


160 000 


80000 


80000 


Ndm 




160 000 


160 000 


160000 


80000 


80000 


Ace. radius [kpc] 


- 


- 


- 


17 


12 


Ace. radial width [kpc] 


- 


- 


- 


10 


5 


Ace. rate 


[Mo/yr] 


- 


- 


- 


5 


5 



gration can result from short-lived density peaks arising from 
interference among density waves overlapping in radius. In the 
near future, the most promising technique to put observational 
constraints on the migration mechanism at work in the IV IW is 
"chemical tagging" (iFreeman and Bland-Hawthomi 120021) with, 
e.g., the HERMES survey, which will allow us to directly mea- 
sure the spread of evaporated open clusters across the Galaxy as 
a function of age (Bland-Hawthorn et al., 2010). 

The effect of bars on galactic disc density profiles 
have been inve s tigated in several previous numerical studies. 
iDebattista et al.l (l2006l) related the occurrence of Type II pro- 
files in barred N-body discs to the bar strength. They have found 
that if the bar is large enough, a break forms roughly at the outer 
Lindblad resonance (OLR), which was inte rpreted as th e result 
of bar-spiral coupling in their simulations. iFovle et all (l2008h 
extended this study further by exploring a more extensive pa- 
rameter space of the dark and baryonic matter in galaxies, with 
combinations of the structural parameters that control the var- 
ious observed profile types and how those profiles evolve over 
time. Both of these studies found that angular momentum redis- 
tribution leads to higher central mass density concentration, to 
the development of a two-component profile and the evolution 
of the inner disc scale length over time toward higher values. In 
both the above works the initial conditions comprised of very ex- 
tended radial disc distributions (such as Type I profiles) in which 
the formation of a break in disc profiles was studied. Thus, they 
could not address the question of what happens when the disc is 
initially limited in extent, as it is most likely the case at the time 
the bar forms. 

In this paper we extend these very interesting results by ex- 
amining the evolution of initially truncated disc profiles as may 
be the case in a more realistic, inside-out disc formation sce- 
nario. Detailed analyses of the disc dynamics resulting in angu- 
lar momentum exchange are presented. We also study the effect 
on disc extension and break formation when initial gas compo- 
nent is present in the disc, as well as continuous gas inflow. 



2. Description of the simulations 

We study three main runs of iso l ated d isc galaxies from the 
GaMer database dPi IVIatteo et al.Ll2007h : the giant SO, Sa, and 
Sbc runs (hereafter gSO, gSa, and gSb). In GalMer, for each 
galaxy type, the initial halo and the optional bulge are modeled 
as Plummer spheres, with characteristic masses Mh and Mb, and 



characteristic radii rn and r^, respectively. Their densities are 
given by 



PhA'') = 



3M, 



H.B 



^K.B 



1 -I- 



-5/2 



(1) 



'H.BJ 



On the other hand, the initial gaseous and stellar discs follow 
Miyamoto-Nagai density profiles with masses Mg and M,, and 
vertical and radial scale-lengths given by bg and ag, and bs and 
fli, respectively: 



P,AR,z)^(^^\. 






(2) 



The parameters for the initial conditions of the three runs, 
including the number of particles used for the gaseous disc, Ng, 
the stellar disc and bulge, A^j, and the dark matter halo, Ndm, are 
given in Table [T] The initial Toomre parameter of both stars and 
gas is taken to be 2 = 1-2 as the initial condition of the Tree- 
SPH simulations an d particle velocit ies are initialized with the 
method described in Hemauisg(ll993l). Fu ll de tails of the simula- 
tions a re given in lDi IVIatteo et al.l ( l2007h and lChihngarian et al.l 
(I2010h . 

Because the gaseous discs are set initially as Miyamoto- 
Nagai density distributions with arbitrary scale-heights, they are 
not in hydrostatic equilibrium. However, equilibrium is reached 
after the first 100-200 Myr We have checked that the time evo- 
lution of the Toomre parameter for gSO (dissipationless) and gSa 
(10% gas) is similar, indicating that the instabilities developed in 
the gSa and gSb models described below are not the result of the 
initial gas instability. 

2.1. Gas and new Stars 

Gas particles are "hybrid particles" (SPH/new stars), character- 
ized by two mass values: one corresponds to the total gravita- 
tional mass of the particle and remains unchanged during the 
entirety of the simulation, while the other is the gas content of 
the particle, decreasing or increasin g according to the local star 
formation rate and mass loss (see IChilingarian et al.l 1201 Ol for 
more details). Gravitational forces are always evaluated using 
the total gravitational mass, while hydro-dynamical quantities 
use the time-varying gas mass. When the gas fraction is below 
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10 
[kpc] 

Fig. 1. Top row: bar morphologies for the gSO, gSa, and gSb models at f = 2 Gyr For each galaxy we show a set of three panels: 
the disc surface density face-on view is presented in the square panel, with the bar side-on and head-on views in the attached right 
and bottom panels, respectively. Bottom row: m = 2 Fourier amplitudes, A2/A0, as a function of radius from f = to f = 3 Gyr 
as indicated by the color bar in the rightmost panel. The formation and evolution of a bar can be seen clearly in all inner discs. 
Deviations from zero in the outer disc indicate the strength of two-armed spirals. 



5% of the particle mass, it is transformed into a genuine star par- 
ticle and the remaining gas mass is distributed over neighboring 
gas particles. We show the gas fraction as a function of radius 
at different times for the gSa and gSb models in Fig. |3] In this 
paper we do not discuss gas and newly formed stars separately, 
therefore, "hybrid particles" and "gas/new stars" will be used in- 
terchangeably. 



2.2. Gas accretion 

In addition to the preassembled discs with an initial gas fraction 
described above, we also simulate the effect of external gas ac- 
cretion. The infalling gas is deposited in an annulus in the disc 
plane, close to the edge of the gaseous disc. At this radius the 
gas is assumed to have already settled in the plane, moving with 
a circular velocity. To cope with the large mass variations of 
the gas component, the number of gas particles is kept constant, 
while the individual particle mass is varied. Unlike the GalMer 
simulation technique described above, where the gas particles 
are hybrid, here the mass which is subtracted from the gas for 
making stars is added to the neighboring disc stars. Therefore, 
in this prescription both the gas and stellar parti cles can increas e 
in mass. For more details we refer the reader to ICombesI (l20 1 1 b . 

We simulate two galaxies with gas accretion: (1) gSaacc, 
which starts with the same parameters as gSa, except for the gas 
being accreted at the constant rate of 5 M^/yi in the annulus 
r = 17 + 5 kpc and (2) gSbacc, same as gSb initially, except for 



the gas being accreted at the constant rate of 5 Mq/jv in the range 
r = 12 ± 2.5 kpc. All parameters can be found in Table [1] 

2.3. The bars in our simulations 

All models discussed in this paper develop central bars. Here we 
show the bars' general properties, such as morphology, strength 
and pattern speed, for the gSO, gSa and gSb models. Further 
analyses and relation to spiral structure are presented in Sec.|6l 

The top row of Fig.[T]shows the bar morphology at f = 2 Gyr 
For each model we present a set of three panels: the disc surface 
density face-on view (square panel), the bar side-on (right) and 
head-on (bottom) views. To make these plots we only considered 
stars in the range shown in the figures, excluding the outer parts 
of the disc. The normalized contour separation levels we use are 
shown in the color bar in upper right panel. A peanut is formed 
in all three cases, quite weak for the weaker gSb bar. It is clear 
that the gSO bar is the strongest of all, followed by the gSa. The 
gSb bar is weaker because the gas fraction is higher, and when 
the gas is driv en to the c enter by t he bar, it actually weakens the 
bar (Bournau d and Com bes. 2002t iBournaud et al.L 120051) . 

The bottom row of Fig. [1] shows the Fourier amplitude, 
Am/Ao, of the density of stars (gSO) and stars + gas (gSa and 
gSb) as a function of radius, where Aq is the axisymmetric com- 
ponent and m is the multiplicity of the pattern; here we only 
show the m - 2 component, A2/A0. Different colors indicate the 
time evolution of A2 radial profile from f = to r = 3 Gyr 
every 150 Myr We estimate these by Fourier-analyzing the 
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0.0 0.5 1.0 



t 



1.5 2.0 
[Gyr] 



2.5 3.0 



Fig. 2. Top panel: circular velocities at f = 2 Gyr for the gSO, 
gSa and gSb models. Middle panel: time evolution of the bar 
pattern speeds, Q^. Bottom panel: m - 2 Fourier amplitudes, 
A2/A0, as a function of time. Note that here we take the maxi- 
mum values inside the bar; these are not sensitive to the lengths 
of the bars which can be seen in the bottom row of Fig. \T\ 



disc surface density d istribution, as described in Sec. 4.11 by 
iDi Matteo et al.l (l2007h . Deviations from zero indicate the pres- 
ence of bi-symmetric structure, i.e., a central bar or 2-armed 
spirals. For example, for the gSO model the bar is identifiable 
by the smooth curve in the inner disc, which at later times has 
a maximum at ~ 3 kpc and drops almost to zero at ~ 8 kpc. 
Deviations from zero seen beyond that radius are due to the spi- 
ral structure. The gSO bar is the strongest for the first half of the 
simulation. At later times its maximum amplitude decrease and 
platoes in the range r = (3,5) kpc, unlike the localized maxi- 
mum at r ~ 3.5 kpc in the case of gSa. The gSb bar is clearly 
much shorter than the other two, peaking at approximately the 
same A2/A0 value as gSa at later times. The time variation of 
the maximum Fourier amplitude values can be seen in Fig. |2] 
bottom panel. At the beginning of the simulation the spirals are 
strong, especially for gSa, decreasing in amplitude at later times 
(blue/black lines) as the disc heats up. At the same time both bars 
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Fig. 3. Gas fraction in the hybrid particles (initially 100% gas) 
at different times for our dissipational simulations. 



extend while their pattern speeds decreases significantly (see 
Fig. [21 This behavior of bar evolution is well kno wn and is likely 
due to angular momentum loss to the inner halo (lAthanassoulal 
l2003h . Note, however, that some of the bar angular momentum 
is also transferred to the outer disc, resulting in the extended disc 
profiles seen in Fig.|4] discussed in the following section. 

The top panel of Fig. |2] presents the circular velocities for 
the gSO, gSa and gSb models in colors and line-styles as indi- 
cated in the figure. We estimate these from the disc stars only 
by correcting for asymmetric drift using the following equation 
(jBinney and Tremaine, 2008): 



Vc 



Vl + crl 



, d(\nv)\ do-^, 

1 H r 

d(lnr)l dr 



(3) 



where V^ is the circular speed and V^ is the measured rotational 
velocity of stars for each bin. The initially identical circular ve- 
locities of gSO and gSa differ at this time inwards of ~ 5 kpc, 
with the gSO starting to decline sooner as the galactic center is 
approached. This is likely due to the larger bar developed in the 
gSO model. Both gSO and gSa have declining circular velocities 
as a result of their large bulges. The gSb circular velocity is flat 
to large radii and similar to what is expected for the MW. 

The middle panel of Fig. |2] shows the time evolution of the 
bar pattern speeds, Q/,, for all models. For the dissipational gSa 
and gSb, Q/, increases quickly after the bar forms and then slows 
down. This results from the mass increase inside the bar due to 
the gas inflow following the bar formation (see Fig.|4] third row). 
The gSb's shorter bar is reflected in its larger patterns speed. 

The bottom panel of Fig.|2]shows the bar strengths as a func- 
tion of time. We get these from the maxima of the m - 2 Fourier 
amplitudes, A2/A0, inside the bar region. Further discussion of 
the bars' pattern speeds, strengths and relation to spiral structure 
can be found in Section|6l 



3. Radial variation of the density and metallicity 
profiles 

In Fig. m we present the temporal evolution of the azimuthally 
averaged stellar disc surface density (first row) and metallicity 
(second row) radial profiles for our galaxy models. The third 
row presents the radial density of the gas and new stars. Bulge 
stars are not shown in order to compare the effect on the initially 
exponential discs for all galaxies. From top to bottom these are 
the gSO, gSa, and gSb models described in Sec.|2] Densities are 
plotted logarithmically, thus an exponential appears as a straight 
line. Different lines in a given panel show each system at f = 
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Fig. 4. Temporal evolution of the azimuthally averaged density (top row) and metallicity (middle row) disc profiles for our galaxy 
models. Bulge particles are not considered. The different times shown are color coded as indicated in the second row, third column. 
The time evolution of the bars' corotation and outer Lindblad resonances are shown by the dotted and solid vertical lines, respec- 
tively, where the colors correspond to the different times. Note that the bars slow down considerably, especially in the simulation 
lacking a gaseous component (gSO). 



0, 40, 60, 80, and 3 Gyr These are color-coded as indicated in 
the third row, rightmost panel. The vertical dashed, black lines 
show the initial scale-length, hd, of each disc. The dotted and 
solid vertical lines give the radial positions of the bar's corotation 
and the OLR for each time indicated by the colors. The time 
evolution of the bar pattern speeds for all models is shown in 
Fig. |2] The initial gas fraction, /gas, is given at the top of each 
panel. 



Examining the density profiles shown in the top row of 
Fig. m we observe that the initial disc particle distributions (red 
solid curves) expand well beyond the initial radius as the discs 
evolve. We note the appearance of two breaks in the initially sin- 
gle exponentials: (1) an inner one at r < l.S/i^ (near the bars' 
CR, dotted vertical lines), inwards of which the profiles steepen 
and (2) an outer one near the bars' OLR (solid vertical lines) 
beyond which the profiles also steepen. The metallicity profiles 
show similar changes with time, except that beyond the outer 
break the gradients reverse. There is, however, a number of re- 
markable differences among the evolution of the three galaxy 
model, which we now discuss. 



3. 1 . The effect of a gaseous component: the gSO and gSa 
models 



Here we compare the gSO and gSa simulations, which have the 
same initial conditions except for the lack of gas in the case of 
gSO. Inspecting the upper left panel of Fig. |4](gS0 model), we 
note that the disc, initially limited to 2.5/ij, expands to over 9lici 
in 3 Gyr There is a break at ~ 13 kpc for earlier times, which 
advances to ~ 17 kpc at the final time, beyond which the ex- 
ponential profile steepens. It is notable that this break does not 
occur at the initial truncation (12 kpc) and, thus, must be related 
to a dynamical process taking place in the disc, rather than the 
initial conditions. At times 0.6 and 3 Gyr, the profile outside the 
bar's CR can be fit by two exponentials reasonably well, with the 
break shifting from ~ 3hd to ~ 4hd. This is consistent with the 
Type II (down-turning) profile, in the classification by EPB08, 
for example. Very interestingly, at f = 3 Gyr (black solid curve) 
the profile between the bar's CR and the outer break appears as 
a single exponential ~ 1 .5 scale-lengths beyond the initial trun- 
cation radius. At this time, the break coincides exactly with the 
bar's OLR, therefore providing a perfect example of the Type 
II.o-OLR profile (see figure 5 by EPB08), where the break in 
the surface-brightness profile happens at an outer ring associated 
with the bar's OLR (e.g., NGC 3504, figure 8 by EPB08). 
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Fig. 5. The distribution of fractional changes in initial radius, (r - r())/ro, for the same time outputs as in Fig. |4] Particles out to 
40 kpc are considered. The distributions are asymmetric across the zero value (vertical dashed line) indicating that the discs expand. 



Inspecting the evolution of the metallicity profile of gSO (first 
row, second column of Fig. HJi, we find that the initially strong 
gradient quickly flattens at r > ly. We observe an upward kink at 
~ 15-20 kpc for early and later times, respectively, similarly to 
the break seen in the density profile but significantly further out 
in the disc. At f = 3 Gyr, the metallicity at r ~ 35 kpc (~ 9/ij) 
is consistent with that at ~ 7-8 kpc of the initial disc (red 
solid curve). This suggests that the stellar population reaching 
35 kpc originates at ~ 7 - 8 kpc, which is near the bar's initial 
OLR and final CR, thus the radius of the strongest migr ation 
taking place in the outer disc dMinchev and Famaevil2010l) . The 
occurrence of breaks in the metallicity profiles is not surprising, 
given the breaks in the surface density and the existence of an 
initial metallicity gradient in our simulations. 

The galaxy model just described lacks a gaseous component. 
The second column of Fig. |4] presents the gSa model, which is 
identical to gSO, except for the 10% gas present initially. In this 
case the breaks in the density profiles at different times span a 
larger range, from ~ 14 to ~ 25 kpc for earlier and later times, 
respectively. Again, very remarkably, the initially truncated disc 
grows as a single exponential between ~ 1 .5/?^ and the break ra- 
dius, this time doubling the disc extent by reaching ~ 6hd before 
it steepens. This model is very much consistent with the Type 
II.o-CT profile, where the break in the surface-brightness pro- 
file is found outside the bar's OLR (e.g., NGC 2950, figure 8 
by EPB08), which has so far been attributed to star-formation 
thresholds. We discuss this further in SecJTTI 

The gSa metallicity profile also exhibits an upward kink at 
later times, but displaced at larger radii, as seen in the density. 

The middle column, third row panel shows the radial den- 
sity profile of the initial gas distribution of the gSa model. At 
later times, these hybrid particles comprise a mixture of gas 
and newly formed stars, which is reflected in their dynam- 
ics. We see the typical behavior expected of gas dynamics in 
barred discs, where bar torques and exchange of angular mo- 
mentum across the bar's CR allow the gas t o flow towards the 
galaxy center (e.g., iFriedIi and Bend 1 9931). yielding the for- 
mation of a "pseudo-bulge" dKormendvl 1 1 9931 : ICourteau et all 
II996l:lKormendv and Kennicutil2004) . Outside the CR, the time 
evolution of the radial density profile of hybrid particles resem- 
bles that of the old stellar component discussed above. A notable 
difference is the break occurring even further out. This is most 
likely related to the larger initial scale-length of the gaseous disc 
(5 instead of 4 kpc). 

How do we explain the differences in the evolution of the 
density and metallicity profiles between the gSO and gSa mod- 
els? Compared to the gSO model, the gSa's bar is similar 



in length but slightly weaker at all times (Fig. |2]l. However, 
the gSa's spirals are stronger, especially during the first Gyr 
(Fig. [U. This is caused by the initial gas component, which 
keeps the disc cool and thu s more unstable. As described by 
iMinchev and FamaevI ( l20IOh . the resonance overlap of bar and 
spirals results in a very efficient exchange of angular momen- 
tum. How effective the migration is depends strongly on the 
strength of both the bar a nd spiral structure (see Fig. 4 in 
iMinchev and FamaevI 1201(31) . Hence, even though the gSa's bar 
is weaker, the presence of stronger spirals results in a substan- 
tially more massive outer disc. 



3.2. The effect of a weaker bar: the gSb model 



The third column of Fig. |4] shows the gSb model, which has 
an initial gaseous disc with mass 20% that of the stellar disc. 
Although the initial stellar disc is truncated at ~ 14 kpc (com- 
pared to 12 kpc for gSO and gSa), we observe only weak disc 
extension in the 3 Gyr time evolution shown in the figure for 
both old stars and gas. The metallicity profile at later times be- 
comes shallower (occurring at ~ 15-17 kpc), as opposed to the 
upturn in the case of the gSO and gSa models. The lack of a 
significant mass transfer outward can be attributed to the much 
weaker bar and spirals. While comparable in strength to the other 
models for a short period of time around 300 Myr, the gSb bar 
decreases in amplitude quickly to a maximum of ~ O.65A2/A0 
(Fig.|2]i and length 4.5 kpc. This is slightly weaker than the gSa 
bar at later times, but significantly shorter. The spiral structure 
seen in the A2/A0 Fourier component (bottom row of Fig. [T} is 
similar to that of the gSO model. Due to the gSb bar's slowing 
down at later times, its resonances move ~ 2.5 kpc outward in 
the disc. This "resonant sweeping" is to be compared to ~ 6 and 
~ 7.5 kpc in the case of gSO and gSa, respectively. As in the case 
of gSa, the gSb density profile evolves to a Type II.o-CT profile, 
but with the break occurring at the initial truncation radius. 

In all three models discussed above, the initial disc scale- 
lengths of the middle exponential (between the inner and outer 
breaks) increase. This is apparent from the flattening of this part 
of the density profiles, where the scale-length is given by the 
negative of the slope of a line fit, since the vertical axis is loga- 
rith mic. This is in agreement with a number of previous works 
(e.g.. lFovle et al.ll2008HD"ebattista et al.ll2006h . 
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3.3. Redistribution of initial radii 

Another way to illustrate the discs' spreading out is by consider- 
ing the redistribution or initial stellar radii. In Fig. Ejwe plot the 
distribution of the fractional changes in initial radius, (r- ro)/ro, 
for the same time outputs as in Fig.|4] Particles out to 40 kpc are 
considered. We see that the distributions are strongly asymmetric 
across the zero value (vertical dashed line), increasingly so for 
later times, indicating that the discs expand. Increase in stellar 
radii can be up to 6 times their initial value for gSO and gSa. On 
the other hand, negative changes in radius are only about 100%. 
The stronger changes for the gSO and gSa models are consistent 
with their more extended discs shown in Fig. |4] 



4. Angular momentum exchange among different 
galactic components 

Before we consider the angular momentum redistribution 
in the disc itself, we first study its exchange between the 
different galactic components: disc, bulge, gas, and dark 
matter halo. Earlier works have related the formation of 
a bar and its strength to the halo concentration ( e.g., 
lAthanassoula and Misi riotis 2002; Bournaud et al.j |2005|) . as 
well as its triaxiality (Machado and Athanassoula, l2010h . As 
bars form, the inner disc looses angular momentum since bars 
live within their CR. Some of this lost angular momentum 
can be acquired by the dark matter Further exchange with the 
halo can slow down bars even more. However, we must note 
that even in simulations wit h rigid halos, bars do slow down 
dCombes and SandersL Il98lh . In that case all lost angular mo- 
mentum goes into the outer disc. Our simulations have live halos 
and we thus expect that some disc angular momentum is lost to 
the halo. 

Following [Di Matteo et al.l(l2008h (their Fig. 4), we show the 
time evolution of the total vertical angular momentum, L~, of the 
different galactic components for the gSO, gSa, and gSb models, 
in Fig.|6] We consider particle distribution within 100 kpc around 
the galactic center. Different colors and line-styles, shown in the 
middle panel, indicate the old stellar component (disc stars), the 
bulge, the gas/new stars, and the dark matter halo. We expect 
exchange of L- among all these components. 

For the gSO model (left panel), the initially non-rotating dark 
matter halo absorbs about 2/3 of the angular momentum lost by 
the disc. Remarkably, the remaining 1/3 is acquired by the (also 



initially non-rotating) bulge. These exchanges are clearly seen to 
take place after the bar is formed and continue until the end of 
the simulation, while decreasing once the bar buckles (and thus 
weakens) at f ~ 1.5 Gyr (see Fig.|2]i. 

The gSa model shows about half the loss in stellar compo- 
nent L- compared to gSO. This is consistent with the larger mass 
transferred into the outer disc we found in the previous Section. 
However, we note that the gain by the halo and bulge (identi- 
cal to those of the gSO model) is very similar to the gSO case. 
The additional angular momentum here comes from the gaseous 
component (dashed green curve), not present in the gSO simula- 
tion. Thus, the gas disc either exchanges L- directly with the non- 
rotating components or passes it on to the stellar disc. Therefore, 
the presence of gas in the disc not only supports stronger spiral 
structure, but also allows the stellar disc to retain most of its an- 
gular momentum, making it available for extending the disc. 

Finally, for the gSb model (right panel of Fig. |6) we see that 
the change in L- for the stellar disc is consistent with zero. Small 
variations are seen for the gaseous disc, the bulge and halo. 
While virtually no disc angular momentum is lost to the halo 
and bulge, the gSb disc was found to extend the least, which is 
related to its weak bar. 

As seen in Fig.|6] only about 10% and 5% of the disc angular 
momentum is lost to the halo and bulge for the gSO and gSa 
models, respectively, and effectively none for gSb. This small 
loss of angular momentum to the hal o is not surprising, given 
that w e have maximum discs (see e.g.. lDebattista and SellwoodI 
Il998h . 

Since in this work we are only concerned with the angular 
momentum in the z-direction, we simply refer to it as L for the 
remainder of the paper 



5. Disc morphology and angular momentum 
exchange 

In Section[3]we have related the extension of discs beyond their 
initial truncation radii and the occurrence of breaks in their den- 
sity profiles to the strength of bars and spirals. We now examine 
the discs' morphology and dynamics in more detail in the effort 
to determine more precisely the causes for the changes we see in 

Fig.m 
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Fig. 7. Columns 1-3: Differential, face-on, density of the gSO, gSa and gSb discs, for three different times. The first three rows 
show the "old" stars, while the last two rows show hybrid particles. All bars are aligned with the vertical axis. Some of the spiral 
discontinuities are labeled by short, solid, red lines. Columns 4-6: Incremental changes in angular momentum, AL, in time intervals 
of 200 Myr, centered on the time of the snapshots shown in the first three columns. Both axes are divided by the estimated rotation 
curve for each galaxy and the total mass in each radial bin, thus the values displayed are approximately equal to galactic radius in 
kpc. The bar's CR, its OLR, and the truncation radius are shown as the dotted-red, solid-blue, and dashed green circles (density) or 
vertical lines (AL) in each plot. The red arrow heads indicate strong increases in AL near the disc break for the gSa model, resulting 
in the more extended disc seen in Fig|4] 



5.1. Discontinuities in spiral arms 

Here we relate the breaks seen in the radial density profiles in 
Fig. |4] to the discs' morphology. In the left three columns of 
Fig.|7]we show density contours for the gSO, gSa, and gSb mod- 
els. The first three rows show the old stellar population at differ- 
ent times, while the last two rows present the gas/new stars of 
the gSa and gSb models. In order to emphasize the asymmetries 
in the outer disc , the quantity plotted i s the d ifferential density, 
£j,/, as done bv lMinchev and OuillenI (l2008h . This is computed 
as (S - l,a.\i)/'^axi, where S is the raw value of the stellar density 
and Sfl ri is the azimuthally averaged density at each time with ra- 
dial bin-size of 0.5 kpc. The times of each snapshot are shown in 



the first row. The bars are aligned with the vertical axis and their 
CR and OLR radii are shown by the dotted red and solid blue 
circles, respectively. Note that resonances can shift inward and 
outward with time if the bar speeds up or slows down. The green 
dashed circles show the location of the breaks seen in Fig. |4] 
For the sake of comparison, for both the stellar and gas discs the 
green circles show the breaks in the stellar density profiles . 

Un like E (not shown here, but see Fig. 1 in iMinchev et al.l 
l20Ilal) . the differential density S^,/ displays strong variations 
with radius with prominent discontinuities (strong decrease in 
amplitude) in the spirals. Some of these minima are shown in the 
Figure by short, solid, red lines bridging the gaps. These gaps in 
the outer spirals appear to move outward as time progresses for 
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the gSa, but not significantly so for the other two models. This 
is very similar to the occurrence of the breaks in the density pro- 
files (green dashed circles in these plots), which we described in 
Sec. [3] Interestingly, the breaks always lie just inside the outer- 
most gaps in the spiral arms, where the spiral pitch angles de- 
crease, as is expected when spirals get weaker 

The discontinuities in the spiral arms outside the bar's OLR 
move outward with time for gSa, but not significantly so for gSO 
and gSb. We showed in Fig. |4]that even at the final time (3 Gyr) 
the gSO disc break occurs at 17 kpc, coinciding exactly with the 
bar's OLR. It, therefore, appears that the slowing down of the 
gSO's bar (which moves the OLR radius outward) produces the 
dominant effect, i.e., this galaxy is bar dominated. Conversely, 
the discontinuities in the case of gSa, which happen at about 
twice the OLR radius at f = 0.8 Gyr, indicate stronger instabil- 
ities in the outer disc. We study in detail the mode coupling in 
the gSa disc in Sec. |7] The strong spirals in this simulation re- 
sult from the increased surface density in the outer regions due 
to the diffusion caused by the bar, while the disc is maintained 
cool by the initial 10% gas present. Although the gSb model 
starts with 20% gas, the disc does not extend significantly. The 
reason for this is the larger initial gas fraction, which prevents 
the formation of a strong bar due to the central mass concentra- 
tion resulting from the gas inflow toward t he cen ter as the bar 
forms. Additionallv. iMinchev and FamaevI (1201 Ol) demonstrated 
that the mixing by resonance overlap of bar and spiral structure 
is a non-linear mechanism, strongly dependent on the strength 
of the perturbers. Therefore, the smaller bar, which forms dur- 
ing the gSb simulation cannot transfer enough mass in the outer 
disc to make it unstable enough to sustain spiral instabilities be- 
yond the initial truncation radius. Another way of looking at this 
is comparing the amount of stellar mass redistributed for each 
model: inspecting the area inclosed between the red and black 
curves, we find a big difference in favor of the gSO and gSa mod- 
els compared to gSb (relatively speaking since the gSb disc has 
a smaller total mass). In turn, due to conservation of angular 
momentum for the entire disc (no significant loss of angular mo- 
mentum to the live halo was found in Section|4|, we expect less 
mass transfer outwards for gSb, and, thus, less disc extension. 

Discontinuities in spiral arms can indicate changes in the 
dominant pattern (transition between inner and outer structure) 
or constructive/destructive interference between different spiral 
modes (different waves propagating at the same radius). Multiple 
patterns are expe cted in a resonance-c oupled system, as recently 
demonstrated bv lOuillen et al.l (1201 ih . where in an N-body sim- 
ulation both the inner and outer spirals waves were found to 
be coupled with the bar For example, the lack of rotational 
symmetry in many of the snapshots shown in Fig. |7]is indica- 
tive of an m - 1 mode, which may be resulting from inter- 
ference between two patterns of different multiplicity (e.g., a 
four- and a three-armed spiral waves). We show that this is in- 
deed the case in Sec.|7] Observationally, lopsidedness has been 
found to be a common p roperty of spiral galaxies. For example, 
IZaritskv and Rixl ( Il997h have shown that ~ 30% of field spiral 
galaxies in a magnitude-hmited sample exhibit significant lop- 
sidedness in the outer discs. Even classical Galaxy analogues, 
such as NGC 891 exhibit big asymmetrie s, clearly resolved a s 
gas accreting along a large filament or arm (iMapelli et al.ll2008h . 

The first three panels of the bottom two rows in Fig. Qshow 
density contours of the gas/new stars component for the gSa and 
gSb models. They appear quite similar to the stellar counterparts 
we just described. The colder population in these plots makes 
it clearer that gaps in the spiral arms indicate the transition to a 
different set of structure (possibly with different pattern speed). 



For example, for the gSa, the arm at {t,x,y) - (0.8, 10, 15) is 
now clearly seen to be independent of the inner nearby one. Two 
additional prominent outer arms as seen at the earlier time as 
well: compare the gaseous and stellar discs at r = 0.6 Gyr The 
case is similar for gSb. 

5.2. Angular momentum exchange at different radii 

We now study the angular mome ntum exchange at differ ent radii 
in our galactic discs. In Fig. 2 bv lMinchev et al.l(l201 lal) we pre- 
sented the changes in angular momentum as a function of radius 
for the duration of the simulations (3 Gyr) forgSO, gSa, and gSb. 
In contrast, in columns 4-6 of Fig. [T] we show the incremental 
changes in the time intervals 0.3-0.5, 0.5-0.7, and 0.7-0.9 Gyr 
We estimate these as AL(r) - L\ (r) - L(r), where L{r) and Li (r) 
are the "initial" (i.e., t=0.3, 0.5, and 0.7) and "final" (i.e., t=0.5, 
0.7, and 0.9) angular momenta as function of radius at each time 
step. Thus, each plot shows how much migration has occurred as 
a function of radius during the time step of 200 Myr, centered on 
the time of the snapshots shown in the first three columns of the 
same figure. To properly display density contours we consider 
only stars with r > 2 kpc. The bar's CR and OLR are shown as 
the dotted red and solid blue vertical lines in each plot. Note that 
in all cases resonances shift outward due to the bar's slowing 
down, not very significantly for the gSb model. The estimated 
angular momentum at each radial bin is divided by the rotation 
curve appropriate for each galaxy model and the total mass at 
that bin. Therefore, the units of L and AL are kpc. This means 
that, for example, the group of stars at (L, AL) x (3,5) is shifted 
radially outwards by ~ 5 kpc. Similarly, stars losing angular mo- 
mentum have their mean (guiding) radii shifted inwards. 

In all three models the regions around the bars' CR radii 
(red dotted vertical lines) display strong exchange of angular 
momentum. The amount of mixing across the CR is approxi- 
mately constant for all times, while the outer disc regions display 
structure which varies with time. As we will show later, this is 
related to the interplay of spiral waves of different multiplicity. 
The amount of migration taking place outside the bar's CR for 
each simulation is consistent with the discs' extension shown in 
Fig. m the gSa disc is affected the most, especially in the outer 
parts. This is reflected in the fast radial advancement of the disc 
break (dashed green vertical line) compared to the other models. 
The red arrow heads indicate strong increases in AL for the gSa 
model, giving rise to the disc extension. Such peaks are not seen 
for the other two models. The constant change in structure in the 
L - AL plane for the gSa disc implies that it is not a single per- 
turber that drives mass to the outer disc but, rather, a network of 
waves. 

It is clear from these plots that to drive the break outwards, 
the disc needs to be unstable in the region near the break. In 
our simulations (not considering gas accretion), this is possible 
if enough mass is transferred outward, while keeping the disc 
cool so that the Toomre instability parameter, Q, ~ cTr/'L, is suf- 
ficiently low. Given the strong bar and spiral structure in the gSa 
model, and the presence of the initial gas fraction, it is natural 
that this disc extends the most. The larger area at positive AL - 
in the gSa outer disc indicates that stars suffer a large increase in 
their home radii (on the order of 3 kpc). The angular momentum 
exchange here is balanced approximately by the smaller scale 
inward migration (~ 1 kpc, red and orange contours) for a larger 
mass fraction. This is in contrast to the gSO and gSb discs, where 
changes above and below the AL = are approximately equal. 

The bottom two rows of Fig. [T] are similar to the above 
plots, but display the gas/new stars component for the gSa and 
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Fig. 8. Changes in angular momentum in the time periods t - 
1 - 2 Gyr (top) and f = 2 - 3 Gyr (bottom) for the gSO, gSa, and 
gSb models. The initial angular momentum L here corresponds 
to f = 1 and t - 2 Gyr, respectively. The time averaged value 
of the bar's CR and 2:1 OLR are shown by the dotted-red and 
solid-blue lines. It is evident that the bars are the most effective 
drivers of radial migration through their CR, despite the fact that 
they are not transient, but only slowly evolving. This is true for 
the entirety of the simulations, as this figure shows. 



gSb models. As in the old stellar population, here we also find 
prominent increases in AL; these are even more strongly empha- 
sized since circular or bits are affected more by the perturbations 
(lMinchevetal.Ll2007h . 

To see whether the bar is an effective driver of radial migra- 
tion for times grater than 1 Gyr, in Fig. |8]we show the changes 
in angular momentum in the time period f = 1 - 2 (top) and 
f = 2 - 3 Gyr (bottom) for the gSO, gSa, and gSb models. 
The initial angular momentum L here corresponds to f = 1 and 
/ = 2 Gyr, respectively. The time averaged value of the bar's 
CR and 2:1 OLR are shown by the dotted-red and solid-blue 
lines. We see that the effect of the bars' CR is evident through- 
out the entire simulations. Even for the gSb model (which hosts 
the smallest bar) changes due to the bar are similar to those of 
an outer spiral between the bar's CR and the OLR. This fact 
renders the bars in our simulation the most effective drivers of 
radial migration, despite the fact that they are not transient, but 
only slowly evolving. This is surprising in view of the theory 
developed by Sellwood and Binnev (2002), where the CR is ef- 
fective onl y if structure i s trans ient. This is consistent with the 
findings of lBrunetti et al.l (1201 Ih . 

To make a better judgement about the drivers of radial mi- 
gration in our discs, we now estimate the pattern speeds. 



6. Pattern speeds 

W e construct power sp ectrograms using the procedure outlined 
bv lOuillen et al.l (1201 ih . Sec. 2.3. During an N-body simulation 
with a live halo (as in our simulations), lopsided modes may 
develop and the bulge and central disc of the galaxy may not re- 
main at a fixed position. We, therefore, subtract the position of 
the centroid of the galaxy bulge prior to computing the Fourier 
coefficients. A sufficient number of time outputs in a given time 
window is required to compute a spectrogram. Since structure 
may be transient and pattern speeds may vary with time, a small 
time-step between outputs is needed for an accurate determina- 
tion of the temporal evolution of density waves in a galactic disc. 



For the simulations presented here we have time outputs every 
10 Myr 

In Fig. |9] we show spectrograms constructed from the m = 
2,4, 1, and 3 (top to bottom) Fourier components for the gSO 
(columns 1-2) and gSa (columns 3-4) models. For each galaxy, 
spectrograms are computed at f = 0.5 and t = 0.8 Gyr with a 
time window Af = 1 Gyr, as indicated in each panel. The ver- 
tical axes show the pattern speed in units of km/s/kpc and the 
horizontal axes plot the galactic radius in kpc. For m = 2 and 
m - 4, the resonance loci are plotted for the CR (solid), 4:1 
LR (dashed) and 2:1 LR (dot-dashed) orange lines, computed 
as mQ, mO. + k/4, and mO. + k/2, respectively. For the m = 1 
and m = 3 cases, we plot niD., mO. + k, and rr£i + k/3. The color 
contours represent the strength (in arbitrary units) of the features 
indicated in the color bar 

Firstly, we note that the waves we see in the spectrograms 
may not relate directly to the morphology of our galactic discs. 
What appeared as well defined spiral arms in the stellar density 
plots of Fig. [Tjmay be structures comprised of superposed waves 
moving at different angular velocities. This will become appar- 
ent later in this Section. 

The structure seen in the m - 2 power spectra (top row of 
Fig. |9|l indicates bi-symmetric waves, such as a central bar and 
two-armed spirals. For example, the gSa bar is identifiable by the 
strongest, fast horizontal feature in the inner disc (between and 
6 kpc). Two additional features (spirals) are present in the outer 
disc at lower pattern speeds. The strongest, slower spiral has its 
2: 1 ILR near the locatio n of the bar's CR, as expected in the cas e 
of mode couphng (e.g., iTaggeret alii 19871 : ISvgnet et al1ll988h . 
By comparing the spectrograms at f = 0.5 and t = 0.8 Gyr, 
we see that both the gSa and gSO bars slows down at the later 
time, in agreement with Fig. |2] While it is easy to discern two 
dominant spiral patterns for both the earlier and later times for 
gSa, structure in the outer disc is not very clear in the later gSO 
plot. This is related to the weak spirals of this model at times 
> 0.5 Gyr, as seen in the bottom panels of Fig. [1] Note that the 
m = 2 spirals in all plots end close to the CR (solid curves). This 
fact, as well as the common decrease in pattern s peed for both 
bar an d spirals is consistent with the findings of lOuillen et al.l 
(I2OIII) . where a dissipationless simulation was studied. The au- 
thors interpreted this as a sign that the system was coupled. 

The m = 4 power spectra are shown in the second row of 
Fig.|9] These represent four-armed spirals, as well as the m - 4 
symmetry in the bar (its first harmonic) at 2Q.i,ar- Due to the fac- 
tor of two higher frequency, the bar pattern speed is much better 
defined in this plot. Given that the bar is mostly bi-symmetric, 
its strength here is roughly half of the m = 2 component. 

The m - 3 Fourier component is shown in the bottom row of 
Fig- ID Two dominant (inner and outer) three-armed spiral waves 
are found to propagate in the gSa disc, in addition to the m = 
2 and 4 structure just described. These are weaker by about a 
factor of 0.5 compared to the amplitudes of the outer, two-armed 
spiral wave, as evident by comparing the corresponding color 
bars. Three-armed patterns in the gSO model are also present, 
but not well defined and weaker, as is the m - 2 structure, due 
to the lack of gas and consequently hotter disc. 

We finally discuss the third panels of Fig. |9] Here we show 
the m - \ Fourier component, which represents lopsided struc- 
ture, such as one-armed spirals. We note that the position of the 
centroid of the inner disc has been subtracted, therefore, the fea- 
tures seen in these plots are most likely due to the interaction of 
different waves at certain radii. It is remarkable that the strongest 
m - \ features in both models (better visible in gSa) coincide 
with the m - 3 spirals in both radial extent and power (see the 
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Fig. 9. Power spectrograms (from top to bottom) of the m - 2,4, 1 and 3 Fourier components for the gSa (columns 1-2) and gSO 
(columns 3-4) models. Spectrograms are computed in a time window Af = 1 Gyr centered on f = 0.5 and t = 0.8 Gyr. For m - 2 
and m - 4, the resonance loci are plotted for the CR (solid), 4: 1 LR (dashed) and 2: 1 LR (dot-dashed) orange lines, computed as fi, 
n + k/4, and D. + k/2, respectively. For the m - 1 and m = 3 cases, we plot Q.,Q. + k, and Q + k/3. Both bars and spirals are seen to 
slow down at later times. The dotter-red and solid-blue vertical lines show the radial location of the bar's CR and OLR. 



color bars). The inner/faster and outer/slower m — \ and m - 3 
waves start at the bar's CR and OLR, respectivel y, hinting that 
these modes may be coupled. [Brunetti et al.l(l201 lb also found an 
m - \ pattern in an N-body barred disc and suggested that this 
feature can be a strong driver of stellar diffusion (migration); the 
nature of this perturbation may be similar to what we find here. 
In the power spectrograms plotted in Fig. |9]we found some 
well defined features, e.g., inner and outer spiral features in the 
gSa disc. Such a low number of coherent structure in spectro- 
grams averaged over a long time period could mean two things: 
(i) spirals are recurrent, but form at the same radii and pattern 
speeds or (ii) spirals are long-lived. In the following section we 
investigate the longevity of spiral structure in our simulations. 



6.1. Spiral structure longevity 

At present the nature of galactic disc spiral structure is not 
well understood. Though it is generally accepted that spirals 
are density waves there exist two competing theories: (i) tran- 
sient/recurrent spirals, and (ii) long-lived, steady-state spirals. 

Recurrent spiral i nstabi lities have been reported by 
ISellwood and Carlberd (1 1 984 and ISellwood andTinl (Il989l) in 
their simulations of self g ravitating discs. It was argued by 
iToomre and Kalnaisl ( 1199 Ih that these transient spiral density 
waves are due to the swing-amplification mechanism as first for- 



mulated bv lToomrel(ll98lh 

Alternatively, the concept of quasi-stationary density waves 
was developed (mostly a nalitically) b y Lin et al. ( 1969) and cul- 
minated in the work by iBertin et aH (Il989a b) and iLowe et aP 
(Il994l) . While thought to always produce short-lived spi- 
rals, N-body simulations have been constructed to yield long- 
lived spiral density waves lasting for over five rotation pe- 
riods, by introd ucing an inner Q-b arrier to shield the 2:1 
ILR (Thomasson et ; 



n g an inner (^-barrier to shield the 1:1 
all 19901: Elmegreen and Thomasson i ll993l : 



iDonner and Thomassoni 119941: IZhangi Il996l) . We note that the 



recent study by ISellwoodI ( 1201 lb have questioned the results of 
these works, claiming that spiral pattern speeds have not been 
measured properly since multiple modes were taken for a single 
one. 

We want to find out how long-lived are the spirals in our 
simulations. As seen in the above introduction, estimating spi- 
rals pattern speeds has been proven to be a hard task. How do 
we ensure we do this properly? Measuring the time evolution of 
the spiral structure seen in the stellar or gas density is not ap- 
propriate since multiple waves (if present, as in our simulations) 
can interfere and thus lead to the wrong result. Power spectro- 
grams, on the other hand, allow for estimating the frequencies of 
the actual waves (at different multiplicities) and is thus, perhaps, 
the most appropriate way. 

In Fig.|9]we plotted Fourier power spectrograms, computed 
using a time window of 600 IVIyr Short transient spiral waves 
cannot be identified in this plot. To look for spiral lifetimes we 
need to decrease the time window used, as well as the time spac- 
ing. Therefore, we construct spectrograms in a time window of 
300 Myr, i.e., half of the one we used before, and plot these ev- 
ery 150 Myr from t - 150 to f = 900 Myr We chose to examine 
the gSa disc because of the better signal to noise ratio, given 
its strong spirals. We present our results in Fig. [10] In contrast 
to Fig. |9] the vertical axes here show the angular frequency, o), 
rather than the pattern speed Q = w/m, in order to facilitate a 
discussion of coupling between diff'erent modes (next Section). 
Because of the additional factor m here, the bar's m - 4 response 
appears at twice the m = 2 w-value (visible at later times as the 
bar slows down). The color bars are normalized for each multi- 
plicity to better follow the variation in amplitudes. 

We first examine the strongest spiral mode, m - 2, shown in 
the first row of Fig. [10] Despite the smaller time window used, 
structure still appears coherent. At f = 300 (first column) a wave 
of frequency u ~ 0.05 Myr"' forms, similar to the outer feature 
we saw in Fig. [9] first column, extending between the bar's CR 
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Fig. 10. Time evolution of the m — 2,4, 1, and 3 (top to bottom) power spectra for the gSa model. The vertical axis shows the 
frequency, o) - niQ. (in units of km s ' kpc '), rather than the pattern speed Q, as done in Fig. |9] in order to facilitate a discussion 
of coupling between different modes. Time outputs are every 150 Myr with a time window Af = 300 Myr Contour levels are 
normalized for each multiplicity to better follow the changes in amplitude. The orange curves show the resonant loci as in Fig.|9] 
It is clear that the outer m = 2 wave has a stable pattern speed (slowly decreasing as the bar slows down) for over 600 Myr, while 
slowly weakening and extending outward with time. In contrast, the inner m = 2 wave bounces between the bar and the outer one at 
the beat frequency of the latter two patterns. The m = 3 waves are also present for the entire time, shifting outwards with time and 
largely coinciding in radius with the m = 1 features. 



and the disc break (see Fig. |4]i. At later times this wave is always 
present, increasing in length up to r ~ 18 kpc at f = 750 Myr 
When animated, this two-armed feature is aways present in the 
spectrograms, smoothly changing from one time output in the 
figure to the next. We, therefore, conclude that this pattern has 
a lifetime > 600 Myr The smooth decrease in strength, as well 
as the radial extent with time, is also in agreement with the con- 
clusion that we see the same pattern in all snapshots. In con- 
trast, the inner m = 2 spiral, while appearing similar to the 
outer one in Fig. |9] we now find it to exhibit strong variations 
with time. When animated, this feature is seen to be driven by 
the interaction of the bar with the outer m = 2 pattern, bounc- 
ing back and forth between the two on a timescale consistent 
with their beat frequency. In other words, every time the (faster) 
bar encounters the (slower) outer spiral wave, this inner wave 
is regenerated, speeding up to catch up with the bar (f = 450 
and t = 600 Myr) and later on slowing down to reconnect with 
the outer spiral (f - 750 Myr). Such an inner structure, con- 
necting the bar with the dominant spiral has been reported be- 
fore and has been proposed to provide an explan ation f or the 



nature of the "long" bar in our Galax y CAthanas soulai 2005 
Martinez- Valnuesta and Ge rhard! l20Ilt iRomero-Gomez et al. 



201 U.Athanassoula,.2012.) . An outer pattern with a roughly con- 



stant pattern speed has al so been reported from simulations o f 
non-barred galactic discs (iRoskar et al.L l20Ill ISellwoodl l20Ilh . 
[Roskar et al. (201 IJ) also saw an inner/faster spiral wave with a 
strong variation in angular velocity (their Fig. 5). 

We estimated above that the lifetime of the two-armed outer 
spiral wave in our gSa model is > 600 Myr. At its average ro- 
tational frequency of w ~ 0.04 Myr \ this corresponds to > 4 
rotations, which is relatively long-lived and comparable to es- 
timates foun d in previous work (e.g ., jThomassonet al. 199' 



Elme green a nd ThomassonlI993tlDonnerand ThomassonlI99 
Zhangil996) . It should be noted that the model we consider here 



has a substantial bar, which may be related to the longevity of 
spirals, especially if mode coupling is present (see next Section). 
When animated, the m = 1,3, and 4 features are seen to 
evolve fairly quickly at the earlier times of the simulation and 
change position and frequency together, sweeping their reso- 
nances through the disc. At later times these waves become 
longer-lived: the m = 1,3 and 4 features at r a: 8 and 16 kpc, 
seen in the fourth column of Fig. [10] form at r a; 50 Myr and 
persist until the final time shown in Fig.[TOl while slowing down 
and moving to larger radii. The coincidence in radial positions of 
all these waves at different times of the disc evolution strongly 
suggest that these modes are coupled. In Section|2]we investigate 
this possibility. 

6.2. Relationship between density waves and spiral arms 

Here we compare the spiral waves found in the power spectro- 
grams to the spiral arms seen in the stellar density plots. We 
turn our attention to Fig. [TT] where the top four rows are the 
same as Fig. [TO] and the fifth row shows differential stellar den- 
si ty contours o f galac tic azimuth versus radius, as was done 



5ity 



b\ IOuilIen et alJ (l20Ilh . Here we see the two halves of the bar 
as the horizontal strong features at low radii and = 90° and 
270°, while spiral structure appears as lines of negative slope. 
Rotation is in the positive vertical direction. The corresponding 
disc face-on contours are shown in the bottom row of the same 
figure. 

An attempt is made to match the waves of different multi- 
plicities to the spiral features seen in the stellar density. Pink 
vertical lines connect the ends of the most likely matches be- 
tween waves and material arms. For most times shown, we can 
identify the slow, two-armed wave and some outer m = 3 and 
m - \ modes. Although for most time outputs the face-on stellar 
structure looks two-armed, we can now see that this is only true 
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for the m - 2 waves up to the discontinuities just outside the 
bar's 2:1 OLR (blue circle). The outer extensions of the mate- 
rial spirals are, in fact, m - 3 and m - 1 modes, easily seen for 
t - 600, 750, and 900 Myr However, a large number of waves 
cannot be identified in the morphology plots because they inter- 
fere. Therefore, one must be careful when making conclusions 
about properties of spiral structure (such as pattern speed and 
longevity) based on morphology alone. 

7. Non-linear coupling 

Multiple patterns in N-body simulatio ns h ave been known to ex- 
ist sinc e the work of lSellwoodI (Il985h and lSellwood and Sparkel 
(Il988l) . who found that a bar can coexist with a spiral pat- 
tern m oving at a much lower angular velocity. iTagger et al.l 
(Il987h and ISvgnet et all ( fl988) explained this as the non- 
linear mode coupUng between the bar and the spiral wave. 
These findings were later con firmed by the numerical stud ies of 
iMasset and Tagged (1 19971) and lRautiainen and Salol(ll999l). 



According to the th eoretical work by iTagger et al.l (Il987b 



and ISvgnet et al.l (Il988l) . two patterns can couple non-linearly 
as they overlap over a radial range, which coincides both with 
the CR of the inner one and the ILR of the outer one. This co- 
incidence of resonances results in efficient exchange of energy 
and angular momentum between the two patterns. The coupling 
between the two patterns generates beat waves (as we describe 
below), also found to have LRs at the interaction radii, resulting 
in a strong non-linear effect even at relatively small amplitudes, 
^autiainen and Salo ( 1999) showed that coupling between a CR 
and 4:1 ILR, as well as ILR and OLR is also possible in N- 
body simulations. Waves couple with a selection of frequencies 
which optimizes the coupling efficiency. Strong exchange of en- 
ergy and angular momentum is then possible among the coupled 
waves. 

A mode with an azimuthal wave number and frequency nii 
and wi, respectively, can couple to another wave with 1112 and clI2 
to produce a third one at a beat frequency. The selection rules 
are 



m — m\ +1712 
and 

U - (jJ\ ± 0)2 



(4) 



(5) 



We next look for relation between the diff^erent waves found in 
our gSa model. 

7.1. Relation between m - 2 and m - 4 



Following iMasset and Taggen (Il997l) . we first search for cou- 
pling between the bar and the slow m = 2 outer spiral wave 
found in our gSa model. Inspecting Fig. [TOJ we find that the 
slower m - 2 feature discussed in Sec. 16.11 has its 2:1 ILR 
coinciding with the bar's CR at all times. This is in agree- 
ment with the expectation that this strong wave is non-linearly 
driven by the bar. In such a case, it is also expected that this 
would result in the generation of two additional waves propa- 
gating at the beat frequency of the partner waves (if they were 
allowed). Using Eqs. |4|and|5]and considering the second col- 
umn of Fig.[TOl we add the corresponding frequencies and wave 



numbers: a)„,4 = w^a,. + a>s 



0.115 + 0.035 = 0.15 and 



Wmo = ciJhar — <^sp ~ 0.08 Myr" . While we do not show here 
m - spectra (ring-like structures in the disc), we can search 
for the predicted m = 4 wave. Such a clump, with oj ^ 0.155 



is identifiable in the m - 4 spectrogram at the same time output 
(t - 450 Myr), in excellent agreement with our estimate above. 
This feature, found with a much better defined pattern speed in 
the longer time-window spectrogram of Fig. |9](left), propagates 
between its 4: 1 ILR and its CR, where it is then strongly attenu- 
ated. Its 4: 1 ILR coincides with the region of resonance overlap 
between the parent waves in accordance with the expectations 
dMasset and Tagg er, 1 997) . 

It is interesting to note that the slow m = 2 spiral wave 
strongly decreases in amplitude near its 2:1 ILR at certain times. 
This is due to the periodic interaction with the inner m-2 spiral 
described in Sec. 16.11 This feature reconnects the slow m-2 
wave and bar every time their phases coincide. Therefore, the 
m - 4 wave regrows on the same timescale. 

7.2. Relation among m - 1,2 and 3 

Previous works have suggested coupling be tween m = 1 
and m-2 modes (Tagger and Athanass oulal 1199 lb . as 
well as among m - 1 ,2 an d 3 modes for both gaseous 
(iLaughlin and Korchaginl [19961) and stellar (iMiller and Smithl 
119921: lOuillen et al.L 1201 lb discs. We now turn our attention to 
the m - 3 spectra, shown in the fourth row of Fig. [TO] For 
most time outputs shown, we see two major clumps extend- 
ing between the 3:1 ILR and 3:1 OLR (orange dashed curves). 
Considering the t - 750 Myr output (fourth row of Fig. [TOt for 
example, we can identify possible coupling among the slower 
m - 2 spiral and the m - 1 and 3 modes centered close to the 
bar's OLR (blue vertical line): (jj„,i = w,„3 - co,„2 ~ 0.07 - 0.04 = 
0.03 Myr '. This is slightly higher than the frequency of the 
m = 1 mode found just outside the bar's OLR at the same time 
output {co,„i ~ 0.02). 



7.3. Relation among m- 1,3 and 4 

With the inner disc's center of mass subtracted, all m = 1 modes 
seen in the third row of Fig. [TOj should be related to coupling 
among patterns with multiplicities m - n and m - n+ \, where 
« is an integer The extent of these features to large radii also 
speaks in support to this claim. By inspecting spectrograms of 
different m-values, we indeed find such evidence. For example, 
considering features overlapping in radius in the spectrograms at 
t - 600 Myr, we can predict the location of m = 1 beat waves 
by subtracting the wave numbers and frequencies of the m - 
4 and m - 3 features at the same time output. For the waves 
centered on the bar's OLR (blue vertical line) and r a; 15 kpc, 
respectively, this results in u),„i - u),„a - Wms ~ 0.095 - 0.075 = 
0.02 Myr-i and w,„i ^ 0.05 - 0.03 = 0.02 Myr"'. These m = 
1 waves are clearly seen at the predicted radial positions and 
frequencies. In addition, a faster m - \ mode is found in the 
same plot very close to the bar's CR moving at w = 0.05 Myr '. 
This one can be explained as the interaction between the weaker 
m - 4 feature with 4:1 ILR near the bar's CR and the strongest 
m - 3 wave. Note the large number of resonance overlaps in the 
bar's CR region: the bar's CR, the m = 4 ILR, the m - 3 ILR, 
and the m-\ CR. 



7.4. Chains of non-linearly coupled waves 

It was already recognized by JLvnden-Bell and KalnaisI (1 19721) 
that trailing spiral structure carries angular momentum outwards. 
Considering a single mode, a wave would absorb angular mo- 
mentum at its ILR and emit it at the CR or the OLR. We may 
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Fig. 11. Rows 1-4 Same as Fig. \TU\ Row 5: stellar density contours of galactic azimuth vs radius. An attempt is made to match the 
waves of different multiplicities to the spiral features seen in the stellar density: pink vertical lines connect the most likely matches. 
Many waves cannot be seen in the stellar disc due to constructive/destructive interference. Row 6: Incremental changes in angular 
momentum computed in 200 Myr spans, centered on the same times as the spectrograms above. The vertical green lines connect 
the likely waves/resonances and their effect on AL. The effect of the inner m - 2 spiral, as it sweeps its CR between the bar and 
the outer m = 2 wave, is shown by the green slanted line in the middle column. Resonant widths increase as the resonance curves 
(orange) flatten in the outer disc: clumps in AL are highlighted by the green shaded rectangles in the last three columns. Row 7: 
Face-on differential stellar density contours for the same times. 



relate this to the large disc extension seen in our gSa model. 
In the spectrograms shown in Fig. [TO]we can identify a continu- 
ously larger number of waves with increasing multiplicity. These 
patterns seem to line up in a specific manner. 

As discussed in Section ItTI the strong, slow m = 2 spiral 
feature has its 2:1 ILR at the bar's CR at all times, suggesting 
non-linear mode coupling between the two waves as described 



bv lMasset and Tagged ( Il997h . In addition to the four-armed beat 
wave resulting from this coupling, on the m - 4 spectrum shown 
in Fig. [To] we also see a major contribution corresponding to the 
first harmonic of the bar at 2(y/,o, (not well shown on this scale, 
since we are more interested in the outer disc structure), as well 
as two other well defined features present at all times. All three 
m - 4 waves with frequencies lower than that of the bar, run 
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Fig. 12. The L-AL plane for a single perturber: a bar (left panel) 
and a four-armed spiral wave (right panel), from the tes t-particle 
simulations presented in iMinchev and FamaevI (l2010h . The bar 
and spiral pattern speeds and strengths are indicated by Q.i,,s and 
eh,s, respectively. The dotted-red vertical lines indicate the CR of 
each perturber and the solid-blue ones - the ILR and OLR. The 
color bar shows the normalized particle density. Such a shape of 
the distribution (indicated by the green line), with an increase in 
AL inside and a decrease outside the CR is characteristic for a 
single perturber, but not when perturbers interact as discussed in 
the text. 



each corresponding column. Due to the large number of modes 
overlapping, it is not an easy task to relate each peak in the 
L - AL plane with a particular resonance region. However, since 
the number of waves decreases in the outer parts of the disc, we 
can associate clumps in the L - AL plane with the m = 3 and 1 
modes seen in the same location (indicated by the green verti- 
cal lines). Because at such distances from the galactic center the 
resonance curves flatten, they run almost parallel to a density 
wave reaching that radius. Consequently, we expect that reso- 
nance radii in the outer disc are much larger than in the inner 
disc. We indeed observe a significantly large, broad clump in the 
changes of angular momentum in the region r = 15-20 kpc, 
indicated by the green-shaded patch. The outermost m = 3 and 1 
density waves seen in the spectrograms have approximately the 
same radial extent as the feature in the changes in angular mo- 
mentum: ~ 5 kpc. Note that this clump becomes broader as it 
reaches larger radii, just as expected if resonance widths were to 
increase. Because we see this effect on a timescale smaller than 
the long dynamical times in the outer disc, we expect that the 
effective resonance here is CR. Remarkably, since the the outer- 
most m = 3 spiral is mostly inside its CR, stars are preferentially 
shifted outward. 



between their 4: 1 ILR and CR, slow down at later times, con- 
sequently reaching larger radii. The 4: 1 ILR of an outer pattern 
is always closer to the CR of an inner one, i.e., there exists a 
chain of CR-4; 1 ILR overlaps, where an inner wave drives an 
outer one by non-linear cou pling, consistent with the findings of 
iRautiainen and Said (1 19991) . 

In contrast to m = 2 and 4, the inner m - 3 mode (and pos- 
sibly the m =1) extend beyond its CR (close to the bar's OLR) 
(although the amplitude does drop outside the CR, especially at 
the later times). Nevertheless, the outer/slower three-armed pat- 
tern has its 3:1 (first order) ILR very close to the CR of the in- 
ner/faster three-armed wave at all times, as seen in Fig. [TO] In 
fact, we observe a third m - 3 pattern in the better defined pat- 
tern speeds at the bottom left panel of Fig. |9](r a; 15 kpc). This 
occurs at the earlier times when the waves are faster, thus allow- 
ing for an additional slower mode. There appears to exist a chain 
of 3:1 OLR-ILR resonance coupling. That figure also allows to 
see better that coupling occurs between the inner waves' CR and 
the outer ones' 3:1 ILR. 

It is hard to search for such coupling in the m - 1 waves due 
to their low frequency and noisier spectrograms. However, they 
are seen to follow the evolution of the m = 3 and m = 4 and most 
likely also interacting. 



8.2. The L-AL plane near the corotation resonance 



We found in the previous section that only the m - 3 wave 
crosses the CR radius, while it still peaks in amplitude inside it. 
The lack of structure extending to the OLR may be relat ed to the 
coupli ng with the bar since it was also found by Ouill en et al.l 
(l2009h . Given this important detail in the type of structure we 
find in our simulations, it is already clear that migration has to 
occur in a different way compared to the mechanism of transient 
spiral structure, where a spiral wave extends between its ILR 
and OLR. As first shown by Sellwood and Binnev (2002), such 
a wave results in a particular shape of the changes in angular mo- 
mentum in the L - AL plane, which is a ridge of negative slope, 
very similar to what is seen near the bar's CR in our plots (which, 
incidentally, d oes not cross its CR). As shown in a number 
of works (iSellwood and Binnevi l2002t iMinchev and FamaevI 
l20I0l:lRoskar et al.Ll20I Ih . this feature has 180° rotational sym- 
metry around the point defined by the intersection of the line 
AL = and the vertical defining the corotation radius. Therefore, 
there is an increase in AL interior to the corotation radius and a 
decrease exterior to it. 



8. Migration mechanisms: the effect of corotation 
versus non-linear mode coupling 

We would like to know the causes for the strong radial migra- 
tion seen in the gSa simulation, for which we find the strongest 
outward transfer of angular momentum. 

8.1. Wide resonances widths in the outer disc 

We showed in the previous Section that the gSa disc presents 
a system of non-linearly coupled modes. As already discussed, 
these are expected to result in strong transfer of energy and an- 
gular momentum at the coupling regions. To see this effect, in 
the sixth row of Fig. [TT] we plot the changes in angular momen- 
tum for the time periods used to construct the spectrograms in 



To m ake this more clear, in Fig. [ 12] we reproduce part of 
Fig. 2 by IMinchev and FamaevI ( |20IO|) and shows the L - AL 
plane for a single perturber: a bar (left panel) and a four-armed 
spiral wave (right panel). The bar and spiral pattern speeds and 
strengths are indicated by Qi_s and et,s, respectively, and corre- 
spond roughly to the ones seen in the gSa model. This shape of 
the distributions, with an increase in AL inside and a decrease 
outside the CR (indicated by the green lines) is characteristic 
for a single perturber If all migration came from the CR, i.e., 
waves were not interacting, then a collection of waves propagat- 
ing at the same time would result in a series of such symmetric 
shapes. However, when non-linear exchange of angular momen- 
tum occurs, such as due to m ode coupling (e.g., i Tagger et aL| 
II987l : lMasset and TaggeJl997l). structure in L-AL may be quite 
different, as shown bv IMinchev and FamaevI (l2OI0l) . 
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Fig. 14. Destruction of the horseshoe orbits near the CR in the vicinity of an ILR or OLR of a secondary wave. First column: 
The red/green particles initially start on rings just inside/outside the CR (dashed black circle) of an m - 2 spiral wave. The top and 
bottom plots are identical. Second column: Same initial conditions, but near the 4:1 OLR (top) or 4:1 ILR (bottom) of an m = 4 
wave. Third column: Same initial conditions, but particles are perturbed by both spiral waves. Note that a secondary wave of only 
1/3 the strength of the first one is enough to disrupt the horseshoe orbits near the CR. Simulations from lMinchev and Ouillenl(l2006h . 



8.3. Changes in angular momentum for cold and hot orbits 

The coupling between waves we found in Section Q and the 
complex structure in the L - AL plane makes it compelling to 
conclude that the migration in our simulations is dominated by 
resonance overlap associated with multiple patterns, rather than 
the pure effect of the CR. However, it may be possible that 
the majority of migration is done by the CR of spiral structure 
but the additional contribution from non-linear mechanisms con- 
taminates the L - AL space. We can check this, since particles 
whi ch migrate due to the CR presumably end up o n cold or- 
bits ('Sel lwood and Binnevll2002l:lRoskar et al.Ll201 ih and, thus, 
should be easy to identify. Therefore, in Fig.[T3]we compare the 
structure in the changes in angular momentum for stellar particle 
samples with different eccentricities at two distinct times for the 
gSa model. Similarly to Figs.lTlandfTTl here in the first row we 
show the changes in AL in 200 Myr centered on 0.6 (left) and 
0.75 Gyr (right). We next estimate the particle eccentricities as 



W 



(6) 



where u and v = Vc - V^ are the residual radial and tangential 
velocities, respectively, while k and Q. are the radial epicyclic 
frequency and the angular frequency. Eq. |6]is the same as Eq. 8 
bylArifvanto and Fuchs (2006), but in slightly different notation. 
For each column, we select subsamples from the total sample 
shown in the first row, for particles with increasing eccentricities, 
e, from zero to 0.5 in units of 0.1, as indicated in each panel. We 
only consider particles with e < 0.5, which comprise ~ 85% 
of the total sample. Therefore, the "total" distribution shown in 
each first-row panel is the sum of all panels below. There are 
no overlaps in the subsets. Note, that we estimate eccentricities 



after the particles have migrated, therefore. Fig. [T3]does not tell 
us whether particles originated from a kinematically hot or cold 
population. 

The panels in the second row of Fig. [13] show particle dis- 
tributions in the range < e < 0.1. As discussed above, for 
these cold orbits we expect to see the well-defined ridges of neg- 
ative slope in the L - AL plane, if migration were dominated 
by the spirals' CR. This expectation is borne out only near the 
bar's CR (shown by the green lines), more symmetric in the 
left plot, although some distortion is apparent at negative AL. 
Nevertheless, for the coldest subsample, the bar's CR is the only 
location resembling what is expected to be the effect of CR. 
This is somewhat surprising in view of the theory developed by 
Sellwood and Binnev. (.2002) . where for effective mixing a den- 
sity wave has to be transient in amplitude so that stars near the 
CR (on horseshoe orbits) remain trapped away from their initial 
guiding radii, once the wave has disappeared. While orbits near 
the bar's CR are similar to those near the CR of a spiral wave, 
our bars remain stable and only slowly evolve with time in both 
amplitude and pattern speed (see Fig. |2]i. As we can see in all 
L - AL plots presented so far, the bar's CR is effective for all of 
our models in the 3 Gyr of the time evolution. This fact renders 
the bars in our simulations the most effective drivers of radial 
migration, despite the fact that they are not transient, but only 
slowly evolving. 

We now examine structure in the L - AL plane for the in- 
creasingly hotter subsamples in Fig. [13] The red, black and green 
arrows in the left column (lined up in the subsequent rows) show 
the location of three of the most prominent clumps with posi- 
tive AL. We find a continuously decreasing fraction of particles 
at the location indicated by the black arrow, with increasing ec- 
centricities, while the opposite effect is seen for the green arrow. 
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Fig. 13. Compaiison between the changes in angular momentum 
for particle samples with different eccentricities at two distinct 
times for the gSa model. Eccentricities are measured after the 
particles have migrated. First row: Incremental changes in AL 
in 200 Myr, centered on 0.6 (left) and 0.75 Gyr (right). Rows 2- 
7: Subsamples from the total sample shown in the first row, for 
particles with increasing eccentricities, as shown in each panel. 
There are no overlaps in the subsets. The green lines of negative 
slope show three CR-like features. The only well defined such 
structure seen outside the bar is found in the fourth row. Note 
that cold orbits in the outer disc (second row) do not exhibit the 
behavior expected at a CR (clearly seen for the bar), indicat- 
ing the effects of non-linear coupling among waves. The three 
prominent peaks shown by the arrows in the upper left panel 
flange differently with variations in eccentricity, suggesting that 
different dynamical processes are responsible for them. 



Interesting, this L-position of the green arrow turns out to be very 
close to the minimum in the distribution at the hottest samples 
(fifth row, left). The peak in AL, indicated by the red arrow, also 
increases with eccentricity and peaks at 0.2 < e < 0.3. Given this 
diff'erent behavior, the three prominent peaks in the top left panel 
of Fig.[T3]are most likely related to distinct dynamical processes. 

An interesting feature is also found at the later time in the 
simulation (right panel of Fig. [T3l l. The (aligned) black arrows 
show the L-position of the densest outer clump seen in the to- 
tal distribution. This feature is progressively shifted to outer 
radii by about 3 kpc as the eccentricities grow form O.I to 
0.5. Additionally, this clump, initially centered slightly above 
AL = 0, moves to positive AL with increasing eccentricities. 
Similar displacement for hot orbits is seen in the first column, as 
well. What this tells us is that particles migrating outward end up 
on hot orbits (compared to the local population). Examining the 
coldest subsamples, we note that the opposite is true for inward 
migrators. This behavior is related to the approximate conser- 
va tion of radial action of migrating stars as discussed in detail 
bv lMinchevet al.l (l20I2ah . In Section l672l we related these outer 
clumps in the L - AL plane to the CR of m = 1 and 3 modes. If 
this is true, then the CR can deliver stars on hot orbits as well, 
contrary to the expectations. 

The only CR-like feature away from the bar's CR is found 
in the left column indicated by the slanted solid green line in 
the fourth row. It is interesting that this structure is smoother for 
orbits with e > 0.1, rather than the cold ones. We now look for 
the wave responsible for that. Inspecting Fig.[TT] we find that the 
nearest CR at this time step (600 Myr) is that of the m - 2 inner 
spiral. As discussed earlier, when animated, this wave bounces 
between the slow, outer m = 2 spiral and the bar, thus sweeping 
its CR through a significant portion of the disc. In the partic- 
ular time depicted here, we happen to catch all its dynamical 
effect: we see it right between the bar and the outer spiral in the 
spectrograms and the time window for making a spectrogram 
is ±150 Myr, i.e., longer than its lifetime. The relation between 
this density wave and the changes in angular momentum at the 
same radius is indicated in Fig. (TT] middle column, by the long 
green vertical line connecting the alleged cause and effect. The 
next nearest CR is that of the inner m - 3 wave, but it is situated 
outside the bar's 2:1 OLR (blue vertical). There is no CR-like 
feature in the changes of angular momentum in Fig. [TT] associ- 
ated with that. We, thus, conclude that the feature outside the 
bar's CR (shown by a green slanted line in both Figs. [TTI and [T3Tl 
is caused by the resonant sweeping of the faster m = 2 wave. 
In other words, neither in the cold nor hot orbits do we find the 
signature of transient spirals of constant pattern speed, which in- 
dicates that interaction of multiple patterns is important, as we 
elaborate further in the next Section. 

While structure in the changes of angular momentum in the 
outer disc varies significantly with changing the particles' ec- 
centricities, the effect of the bar's CR seems to be similar in all 
plots, i.e., the bar's CR is able to migrate stars on both cold and 
hot orbits with similar efficiency. 

8.4. Destruction of the horseshoe orbits near the corotation 
resonance 

We showed in Fig. [12] that near the CR both the bar and spiral 
waves (as single perturbers in test-particle simulations) give rise 
to a line of negative slope in the L - AL plane. However, we 
showed in Fig. [13] that even for cold orbits this shape was not 
seen for any of the spiral waves, except possibly for the one with 
the variable pattern speed (discussed at the end of the previous 
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Fig. 15. First row: Time evolution of the radial velocity dispersion profiles, cr,.(r), for the gSO, gSa, and gSb models. Different 
colors and line styles correspond to the times indicated in the second row, first panel. The bars' CR and 2:1 OLR are shown by 
the dotted and solid vertical lines, respectively. Second row: Same but for the vertical velocity dispersion, crr{r). Third row: Same 
but for the ratio cr~/o-,.. The slope in cr-fcr,. outside the bar shows that the radial profile flattens more with time than the vertical 
one, while its overall decrease means that the disc heats more radially than vertically. This is related to the conservation of vertical 
action. See text for discussion. 



Section). What could be the reason for the lack of this expected 
behavior? It is easy to realize that the shape of structure in the 
L - AL plane is related to the shape of the particles' trajectories. 
Therefore, we now look at the behavior of these "horseshoe" 
orbits near the CR in the presence of only one and then two spiral 
waves. 

In a study of the efi'e ct of multiple spiral waves on the dy- 
namics of galactic discs, iMinchev and OuillenI (l2006h showed 
that in spatial regions where resonances overlap, particle mo- 
tions become stochastic. This is illustrate d in Fig. \T4\ where we 
reprod uce a combination of Figs. 9-11 bv iMinchev and OuillenI 
(|2006|) . The strength of the spiral waves i s given by the values 
of e in the plots, in units as described by 'Minchev and Ouillei^ 
(|2006|) . Both of these are weaker than the estimated strength of 
theM W spiral structure (see discussion in IMinchev and Famaevl 
I2010h . 

The first column shows a face-on view of particles starting 
initially on a ring just inside (red) and a ring just outside (green) 
the CR (dashed black circle) of an m - 2 spiral wave of an in- 
termediate strength. Top and bottom plots are identical for the 
first column. The efficiency of the CR to migrate stars is im- 
pressive: stars from the inner (outer) ring are shifted outwards 
(inwards) by about 3 kpc. Note the unbroken, banana, or "horse- 
shoe", shape of the orbits. 

The second column presents the same initial conditions, but 
near the 4:1 OLR (top) or the 4:1 ILR (bottom) of an m = 4 
wave. Here only small distortion is apparent, due to (1) the week 



spiral amplitude and (2) the fact that near Lindblad resonances 
stars mostly increase their eccentricities by hardly aff'ect angular 
momenta. 

Finally, in the third column the two rings of particles are per- 
turbed by both spiral waves shown in each row. Both top and 
bottom panels are in the reference frame of the corotating, m - 2 
spiral. Remarkably, the secondary wave, which has an amplitude 
of only 1/3 that of the first one, is enough to disrupt the horse- 
shoe orbits near the CR. Moreover, the distances travelled by the 
particles are now much larger, especially in the first row. The 
irregular shapes of these orbits suggests the presence of stochas- 
ticity. Such disruption in the horseshoe orbits can naturally ex- 
plain the lack of structure in the L - AL plane, expected for the 
pure effect of CR. 

We have just shown that a secondary, weaker spiral, is al- 
ready enough to increase dramatically the effect of CR of the 
primary, corotating spiral. In our simulations we have not two, 
but a much larger number of waves overlapping at any given 
radial range. We conclude, therefore, that the strong migration 
seen in our simulation comes from non-linear effects, rather than 
the sole eff'ect of the CR. 



9. Velocity dispersions and disc thickness 

Previous theoretical studies that focussed on radial migration 
have not discussed the evolution of disc velocity dispersion pro- 
files nor disc thickness. However, knowledge of these properties 
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is important as the spiral and bar instabilities, which drive the 
mass build-up in the outer disc, are also expected to be effective 
at disc heating in both radial and vertical directions. When per- 
turbations at two (or more) pattern speeds, such as a bar and a 
spiral density wave are present in the disc, the stellar dynamics 
can be stochastic, particularly near resonances associated wit h 
one of the patterns dOuillenl l2003t iMinchev and OuillenL l200d) . 

In Fig. \T5\ we show the time evolution of the radial, cr^(r), 
and vertical, cr-(r), velocity dispersion profiles for the old stellar 
populations of the gSO, gSa, and gSb models in the first, second 
and third columns, respectively. Initial bulge components are not 
considered. In each panel the different color curves show differ- 
ent times, as indicated in the bottom left panel. The dotted and 
solid vertical lines in the first row show the radial positions of 
the bar's CR and OLR with colors matching the different times. 

We first compare the cTr profiles of the gSO and gSa models, 
reminding the reader that these two simulations start with the 
same initial conditions except for the initial 10% gaseous disc 
for gSa. We note some interesting differences. Firstly, the gSa's 
disc is hotter at all radii and all times after bar formation, sig- 
nificantly so in the inner disc. The stronger heating inside the 
bar's CR of gSa is due to the central mass concentration com- 
ing from gas inflow ins ide the bar following its formation, (e.g., 
iFriedIi and BenzlI993L see also the gas density profile, third row 
of Fig.|4]i. While this simulation has a weaker bar than the dis- 
sipationless gSO, its stronger spirals (see Fig. |7) results in the 
larger increase of cr^ in the outer disc, as well. It is notable that 
the vertical velocity dispersions increase much less than the ra- 
dial ones (true for all models). For instance, contrast this differ- 
ence for gSOatr~ 12 kpc: ~ 45 km/s (a factor of three) increase 
in (Tr and ~ 12 km/s (a factor of two) increase in cr-. 

A remarkable feature in both the gSO and gSa radial veloc- 
ity dispersion profiles is the flattening at later times with pivot 
points at their CR. Similar changes are also seen in the vertical 
velocity dispersion profiles, however the effect is significantly 
stronger for crr{r). This is evident in the ratio cr-,/cr,. shown in 
the bottom row of Fig. [15] where a decline with radius is seen, 
especially for gSO. We have recently described this phenomenon 
injMinchev et al. (2011b, 2012a), relating it to the conservation 
of radial and vertical actions of particles as they migrate. In a 
simple galactic disc model these can be written as J~ - E,/v 
and Jp = EpJK, respectively. As stars migrate, they roughly pre- 
serve their actions but not their corresponding energies, which 
vary with radius as the radial and vertical epicyclic frequencies, 
K and V, respectively. Thus, as stars migrate inwards (outwards) 
they heat up (cool down). Stellar radial migration also mixes the 
radial distributions of the actions {Jp) and {J-): the more the disc 
mixes, the weaker their variation with radius. Therefore, even if 
the migration rate remains constant, for a coeval stellar popula- 
tion we can expect less heating due to migration at later times. 
We can see that this is indeed the case by observing that in the 
region of strongest migration, r < 2.5hd ~ 10 kpc, velocity dis- 
persion increase saturates quickly with time, while the disc out- 
skirts continuously heat since all the stars populating them are 
outward migrators (with larger actions). 

The gSb model also exhibits flattening in the cr,. and cr, pro- 
files, but not as pronounced due to its weaker bar and spirals and, 
thus, mixing. 

We now consider the disc thickening in the three models 
discussed so far. In Fig. [16] we display edge-on density con- 
tour plots for t = 0.4 and 3 Gyr in the left and right columns, 
respectively. The CR and OLR are indicated by the dotted red 
and solid blue vertical lines. Note that as the disc extends it also 
thickens considerably, as expected from the increase in cr, we 
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Fig. 16. Edge-on stellar density contours for the gSO, gSa, and 
gSb models at f = 0.4 and 3 Gyr The dotted-red and solid-blue 
vertical lines show the bars' CR and OLR. Note that as the discs 
extend they also thicken. 



showed in Fig. [15] The increase in scale-height at 2-3 disc scale- 
lengths is double the initial thickness (~ 300 pc) for gSO and 
gSa and increases by a factor of 1.5 for gSb. Note that the disc 
subjected to the strongest migration thickens the most (gSa in 
this case) as expected from our earlier discussion in this Section. 
Nevertheless, the disc thickness is not sufficient to be considered 
a thick disc, with the additional inconsistency of flaring at large 
radii. This has led Minchev et al. ( 2012a) to the conclusion that 
radial migration cannot resul t in the formati o n of thi ck discs and 
most recently confirmed by iMinchev et al.l (l20I2bl) for a MW- 
like simulation in the cosmological context. 

An important consequence of the results of this Section is 
that if extended discs originate through migration, then disc out- 
skirts must be kinematically hot. 



10. The effect of gas accretion on the radial disc 
profiles 

So far in this paper we have neglected the effect of mergers and 
gas accretion by filaments. Studying the effect of these would 
require to explore a large variety of galaxy formation models 
in a cosmological context, investigating how the results would 
depend on changing the rate and radial extent of accretion. While 
we defer such a detailed study to another work, in this Section we 
briefly explore how our results change when smooth, in-plane 
gas accretion is considered for the gSa and gSb models, which 
we refer to as gSa^cc and gSb(,f.c, respectively. Description of the 
simulation technique was given in Sec. 12.21 and parameter for 
both runs can be found in Table [T] We present our results in 
Fig. [17] Starting at the initial time, gas is being accreted at the 
constant rate of 5 M^/yi in the range indicated by the pink strip 
in each column. Varying the extent of this radial range does not 
change the qualitative results. 

Fig. [represents the changes in the azimuthally averaged stel- 
lar density profiles. Initially the discs extend in Type II profiles, 
similarly to the gSa and gSb models, shown in Fig. [4] However, 
after ~ 2 Gyr of gas inflow, the outer stellar discs develop 
Type III (up-turning) profiles. While the stellar density shown 
in Fig. [T7] contains some contribution from the newly formed 
stars, we have checked that the effect on the old stellar popula- 
tion alone is negligible. We also note that the total (stars + gas) 
density profiles (not shown) are very close to Type I. The stars 
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Fig. 17. The effect of smooth, in-plane gas accretion. The pink 
strips show the gas accretion region for simulations starting with 
the initial conditions used for the gSa and gSb models. The ra- 
dial disc profiles are shown in the first row. An up-turn in the 
radial velocity dispersion profile, cri(r), close to the disc break is 
evident in the second row, while the vertical velocity dispersion 
continuously decrease with r (bottom row). 



and gas appear to rearrange in such a way as to preserve a single 
exponential for the total baryonic density. 

The middle and bottom rows of Fig.[T7]show the effect on the 
radial, crr(r), and vertical, o-^{r), velocity dispersion radial pro- 
files. It is remarkable that crr{r) develops an up-turn in the vicin- 
ity of the outermost break in the density profiles. In contrast, cr^ 
at all times continues to decrease with radius. This indicates that 
the stars forming the disc outskirts are on highly eccentric orbits, 
but well confined to the disc plane. We have to note that this may 
be a result of the in-plane accreting we consider here. 

11. Comparison to observed radial disc profiles 

11.1. Type I profiles 

Type I profiles are single-exponential profile. In the classifica- 
tion by EPB08, these extend to at least four disc scale-lengths. 
In terms of bar length s, the statisti c s for /"break /^bai-(^ a.mn.v) for 
the SO -Sb galaxies of lErwin et all (|2008|) and iGutierrez et al.l 
(1201 Ih . are Median = 4.6 and Mean = 5.7 ± 2.5 (min = 3.3, 
max = 14.4). In our simulations a single exponential can be fit 
to the gSa radial profile out to ~ 25 kpc (Fig.|4] top row, middle 
column), which is about five disc scale-lengths at the final time, 
or ~ 5.5 bar lengths, considering the gSa Ao^max ~ 4.5 (Fig.[T]). 
Therefore, our gSa model can be considered a case of Type I. 

However, we note that while it is possible that the gSa profile 
could be mistaken for a Type I in the limit of a low S/N image 
that prevents measuring the profile very far out, there still remain 
a number of Type I profiles extending to much larger radii, which 
cannot be explained with this model. 

It has been pointed out bv lSchava (l2004t) that the main prob- 
lem for the star formation threshold as an explanation for profile 
breaks is that it fails to explain purely exponential galaxies out 



to large radii. While here we show that radial mixing can be the 
driver for these, it should be kept in mind that stars arriving at 
the disc outskirts by migrating over large distances will be kine- 
matically hot as shown in Fig.fTSJ 

11.2. Type II profiles 

All of the galaxy models presented in Fig. |4] are consistent with 
Type II profiles: a single break beyond which the exponential 
steepens. There are three observed flavors of this type, as de- 
scribed by EPB08: 

(i) The break can happen at the bar's OLR (Type II.o-OLR). An 
example of this is fou nd in our gSO model in agreement with 
iDebattista et all (l2006h . 

(ii) The disc break has also been observed to occur outside the 
bar's OLR (Type II.o-CT). These profiles a re currently only 
explained by star-formation thresholds (e .g., lKennicutgll989t 
ISchavell2004t lElmegreen and Hunte3l2006l) . We now show that 
our gSa model results in such a profile naturally, as a result of 
the strong angular momentum outward transfer It is impressive 
that the break is more than two scale-lengths beyond the bar's 
OLR after three Gyr of evolution. The final state of the gSb disc 
also exhibits a Type II.o-CT profile, occurring at the initial disc 
truncation radius. 

(iii) The third kind is the Type II. i profile, for which the break 
occurs at or interior to the end of the bar We do not find such an 
example in our simulations. 

11.3. Type III profiles 

Type III profiles exhibit a shallower exponential beyond the 
break radius, i.e., the disc density drops-off slower with radius 
than in the inner disc (the disc is "anti truncated"). One possibil- 
ity is that these result from gas-rich, minor mergers, as shown by 
lYoungeretaII(l2007h . 

We now show for the first time, that Type III profiles can also 
form in a smooth gas-accretion event, such as from cosmic fila- 
ments (Fig.fTTli. The gSaacc appears consistent with a pure Type 
III, although both the inner and outer segments are not fitted very 
well by single exponentials. On the other hand, gSb„fc is a case 
of a Type II + III, composite profile. However, it is clear that, 
depending on the gas-accretion radius and properties of the bar 
and spiral structure, a pure Type III profile can result. We, there- 
fore, suggest that observed Type III profiles, combined with high 
velocity dispersion in the up-turning region, could be a signature 
of gas accretion. 

An interesting result from our gas-accretion models is the 
Type III and Type II breaks for the existing stellar population and 
the gas, respectively, at the accretion region, while the overall 
baryonic density remains roughly exponential. It is possible that 
this can later evolve into a Type I or a Type II stellar profile. 

The correlation between Type Ills and unbarred or weakly 
barred galaxies (EPB08) is consistent with the model pre- 
sented here, as it is expected that bars are weakened 
e.g., citealtathanassoula05,debattista06) and even destroyed 
iBournaud and Combesll20()2l : ICombesll201 ll) in a gas accretion 
event. 



11.4. Composite profiles 

Composite profiles are profiles composed of several different 
types. These are not uncommon. For instance, EPB08 and 
IGutierrez et al.l (1201 ll) found 6% (barred) and 8% (both barred 
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and unbarred) of there samples, respectively, to exhibit Type 
II+III p rofiles, e.g., NGC 398 2 in EPB08, or NGC 3455 and 
5273 in [Gutierrez et al.l (l2011b . An example of a Type II+III in 
our simulations can be found in the gSba„. model (Fig.fTTI). We 
again stress the fact the stars beyond the outermost disc break 
are expected to be kinematically hot (see Fig.[T7l third panel). 

12. Discussion and conclusions 

We have investigated the time evolution of initially truncated 
galactic discs, via Tree-SPH N-body simulations. We find that 
due to radial migration and torques associated with bar and spi- 
ral instabilities, discs can triple their initial extent, pointing out 
to a strong annular moment outward transfer. If the bar and spiral 
structure are sufficiently strong, a single exponential can extend 
up to ~ 6 initial scale-lengths (~ 4 bar lengths), i.e., doubling 
the initial disc radius, before it steepens. If resulting from radial 
migration (as opposed to an in-situ formation), old stellar pop- 
ulations in disc outskirts must exhibit high velocity dispersions, 
especially in the radial direction, and thus a thickened disc com- 
ponent. Observations of edge-on galaxies, or measurements of 
the vertical velocities of face-on systems, can be used to distin- 
guish between these two different possibilities. 

In the presence of a central bar, the transfer of stellar mass 
outward in the disc in the first two disc scale-lengths is domi- 
nated by the diffusion resulting from the bar's CR and OLR inter- 
acting with resonances associated with the inner spirals (always 
present in a barred disc). We have shown that the bar's corota- 
tion is the most effective driver of radial migration for the extent 
of our simulations (3 Gyr). For additional migration outwards, a 
chain of spiral waves, reaching beyond ~ 3lid is needed for an- 
gular momentum transport. Self-gravitating spirals at such large 
radii can result by increasing the stellar density in that region, 
provided the bar is strong enough (or by gas accretion). A strong 
bar alone, however, is not sufficient. The stellar population in 
the outer region needs to be cool enough so that the Toomre in- 
stability parameter, Q, ~ cr,.IY,, is sufficiently low. Indeed, even 
though the gSO model develops stronger bar, the lack of gas re- 
sults in large velocity dispersion in the region where the outer 
spirals are needed to carry the angular momentum further out . In 
contrast, the 10% gas fraction present initially in the gSa model 
is enough to prevent this from happening, and thus, spirals ex- 
tend to radii greater than 4/1,/. Chains of non-linearly coupled 
spiral waves can develop and carry angular momentum much 
further out than a single one could (see Section|7]i. Note that the 
gaseous discs are not initially in equilibrium, but this is reached 
after 100-200 Myr. Even if the disc is initially not completely 
stable, it is not this initial instability that determines the evo- 
lution and migrati on we discuss in the pa per. An indication of 
this is contained in lMinchev et al.l (l20IIah . where we examined 
how the disc dynamics change with increasing the resolution. 
We showed that although asymmetries developed later (well af- 
ter the initial relaxation) when a higher number of particles was 
employed, the overall results were in agreement. This means that 
independently of the initial relaxation of the system, it is the ap- 
pearance and presence of strong stellar asymmetries that drives 
migration (also in agreement with Brunetti et al. 201 1). 

In the simulations without gas accretion, the single initial 
exponential disc profile quickly turns into three segments fol- 
lowing the formation of a central bar. These can be described 
roughly by three exponentials (Note that all plots in this paper 
exclude the initial bulge component). 

(1) The first exponential is formed by steepening the density 
within the disc scale-length, hd- Increasing the density in the 



first scale-length is equivalent to building a bulge-like compo- 
nent and is related to the buckling of the bar. The velocity disper- 
sion in those regions are also consistent with bulge kinematics 
(Fig.fTSl). This mass transfer toward the central region is caused 
by bar torques, angular momentum exchange across the bar's CR 
and interaction of spiral and bar resonances. 

(2) The second segment forms by flattening the initial density 
profile and ends either at the initial truncation radius (gSb), the 
bar's OLR (gSO), or well outside the bar's OLR (gSa). This ex- 
ponential continuously flattens with time, i.e., the scale-length 
increases, as this is where the strongest disc asymmetries occur 
resulting in the strongest effects of radial migration. We related 
the outer break of this exponential to discontinuities in the spiral 
structure and to significant changes in angular momentum due 
to multiple patterns and torques (Sections |5]and|7]i. 

(3) The third exponential starts either at the initial truncation ra- 
dius (gSb) or further out (gSO and gSa), as described above. It 
can extend to ~ lO/j^ if migration is strong enough. We find 
that the slope of this disc segment is a strong function of time, 
as is the location of the break at which it starts. The longer the 
disc is subjected to spiral and bar perturbations, the less steep 
it becomes (all our models) and/or the farther out the break is 
displaced (gSa model). A notable feature in outer discs resulting 
from radial migration is that stellar populations are hot, particu- 
larly in the radial direction. 

We note that a fourth exponential can result when gas accre- 
tion is considered. This outermost segment is shallower than the 
inner two, as in a Type III profile (see Fig.fTTTi. 

To find out what gives rise to the strong outward transfer 
of angular momentum we examined in detail the disc dynam- 
ics. We first showed that outer breaks and discontinuities in 
the spiral structure seen in the density plots coincide with sig- 
nificant changes in angular momentum occurring at the same 
radii (Section |5]i, suggesting the interplay of resonances associ- 
ated with multiple patterns. By creating power spectrograms, we 
showed that multiple waves of multiplicities m - 1,2,3 and 4 
are indeed present in our discs (Section |7), with strong evidence 
for non-linear coupling, as seen in a number of previous works 
(Tagger et al., 1987; Svgnetetal., 1988; Rautiainen and Sal3, 
ll999tlMasset and Taggeiill997HOuillen et al.Ll20Ilh . Thus, the 
various waves involved conspire to carry the energy and angu- 
lar momentum, extracted by the first mode from the inner parts 
of the disc, much farther out than it alone could. In contrast 
to strong turbulence, where a large number of modes interact, 
mode coupling occurs in situations where only a small number 
of waves or modes can exis t, so that each in teracts non-linearly 
with only a few others (e.g.. lDavidsonlI972l) . This small number 
of active modes translates into long correlation times, i.e., the 
quasi-stationary structure found in the simulations (> 600 Myr). 

At radii larger than ~ 15 kpc the curves defining the reso- 
nances in the power spectrograms (Figs. 1911111 become increas- 
ingly flattened. Therefore, the resonance widths there must be 
larger than in the inner disc. This may have profound conse- 
quences on the formation of extended disc profiles, as stars will 
be affected over a large range of radii. For these outer regions, the 
outer spiral's CR of our gSa model spans a wide range of radii, 
giving rise to effective migration. Remarkably, since this spiral 
wave is mostly inside its CR, stars are preferentially shifted out- 
wards, as is expected to happen inside the CR. We note that the 
breaks in the density profiles in our simulations app ear near the 
CR of spirals in agreement with'Roskar et al.l(l20I ij). r ather than 
the spirals' OLR, as suggested by Debattis ta et al.l (12006). 

We showed in Section[8]that the incremental changes in an- 
gular momentum (the L - AL plane) exhibit complex behavior. 
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unlike what is expected from the effect of the corotation reso- 
nance (CR) alone. Transient spirals have been suggested to mi- 
grate stars without increasing their eccentricities. In an effort 
to separate their effect from that of resonance overlap, we di- 
vided the gSa disc into different stellar subsets according to their 
eccentricities. We expected to find well defined ridges of neg- 
ative slope in the L - AL plane for the cold population, were 
the disc mixing dominated by the CR of spiral waves. However, 
we found no such features but only the complex, clumpy struc- 
ture related to the non-linear coupling among the bar and spiral 
waves of different multiplicitie s. We explained this w i th the dis- 
ruption of the horseshoe orbits (iMinchev and Ouillenl l2006b ex- 
pected near the CR, from resonances associated with other den- 
sity waves. We conclude that although the CR is very effective 
at mixing, in the presence of more that one density wave in the 
disc, one cannot expect a linear superposition of the correspond- 
ing CRs, but rather, a network of waves working together 

The effect of the bars' CR in our models is evident through- 
out the entire simulations time of 3 Gyr (Figs.|7]and|8). This fact 
renders the bars in our simulation the most effective drivers of 
radial migration, despite the fact that they are not transient, but 
only slowly evolving. This is surprising i n view of the theory 
developed bv lSellwood and Binnev (! 2002h . where the CR is ef- 
fective only if structu re i s tra nsient, but agrees w ith the work by 
iBrunetti et ail (l2011h and lMinchev et"aLl j2012bl) . 

We also considered the effect of gas accretion (Section [2. 2b 
and found that this results in the formation of Type III profiles 
in the stellar population. The stars and gas appear to rearrange 
in such a way as to preserve a single exponential for the total 
baryonic density. 

iFovle et al] (l2008h found that once a two-component profile 
is established, (i) the break remains at the same physical loca- 
tion over time and (ii) the outer profile preserves the initial dis- 
tribution slope. The time for break development depends on the 
ratio nidi A, where m^ is the ratio of the disc mass fraction and A 
the halo spin parameter In contrast to (i) above, in the present 
work we find strong temporal dependence for our gSa and gSO 
models, where the break advances considerably with time over a 
period of 3 Gyr We can reconcile this difference by notin g that 
our discs evolve only to 3 Gyr, while iFovle et all (l2008b study 
their galaxies to 10 Gyr We emphasize that the spiral structure 
is very important for the outward transfer of angular momen- 
tum. In the absence of gas accretion, spiral instabilities weaken 
considerably in the outer discs due to the increase in velocity dis- 
persion, thus preventing further evolution in the disc profiles. For 
example, while Foyle et al. did not find a break in their model 98 
(their Fig. 17) although a bar was formed at f ^ 6 Gyr, the reason 
for this is most likely the lack of sufficiently strong spiral struc- 
ture (at that later evolution time) needed to transport stars to the 
outer disc. Similarly to our findings, for their model 41 (Fig . 16) 
the break can be seen to evolve from r » 2.5hd at t — \ Gyr to 
r ^ A. 5hd at r = 4 Gyr Another difference between [Fo vie et all 
(l2008h and our work is that we find that the outer disc profiles be- 
come flatter with time (ii above). This discrepancy is most likely 
related to the lFovle et al.l (l2008h 's initial density distributions ex- 
tending to larger radii, in contrast to our initially truncated discs. 
Therefore, changes in their preexisting outer disc profiles may 
be hard to see. 

Our findings are c onsistent with iPerezI (l2004b and 
iTruiillo and PohlenI (l2005h . who have shown that observed face- 
on galaxies at high redshift (z ~ 1) have breaks at smaller galac- 
tocentric radii relative to nearby galaxies. Of course, this can 
also be related to the disc inside-out growth. 



Ivan den BoschI (1200 ih assumed that the distribution of ob- 
served (inner) disc scale-lengths would be coupled to the spe- 
cific angular momentu m of the halo. It was pointed out by 
iDebattista et alJ (l2006h that the inner disc scale-length evolves 
considerably over time, as a result of star formation, and cannot 
be assumed to be representative of the initia l profi le disc scale- 
length. While agreeing with this. iFo vie et al.l(l2008h have argued 
that the outer disc scale-length remains stable over the discs evo- 
lution and is a close match to the initial disc scale-length. In view 
of our findings, where strong changes with time are observed in 
both the inner and outer profiles, we conclude that both the inner 
and outer scale-lengths cannot be used as good tracer of the spe- 
cific angular momentum of the halo. Note that this conclusion 
is based on our simulations performed in isolation. It remains to 
be seen whether this is still a true statement in the cosmological 
context, where there would be a continuous supply of new gas 
accreting all over the disc. 

Stars migrate along spiral arms and are deposited in the disc 
outskirts on a timescale consistent with the pattern speed of the 
outer spiral density wave. This process will progressively deposit 
increasingly more metal-rich stars from the inner disc in the out- 
skirts. Given the long dynamical times at large distances from 
the galactic center, we expect that strong azimuthal variations in 
metallicity might exist in the outer parts of galactic discs. This is 
the subject of a follow-up work (Minchev et al.. In preparation). 

While we investigated the effects of secular evolution in 
barred discs and the effects of gas accretion, our models con- 
sist of pre-assembled discs. For example, a mass increase due 
an accreting thin disc would shrink the older, thicker compo- 
nent and decrease the velocity dispersion we find in the disc out- 
skirts. Therefore, our results may not accurately represent ob- 
served galaxies. Rather, our study will be useful in separating 
the effects of pure dynamical evolution from those of cosmo- 
logical environment and merging histories. The tidal effects of 
small s atellites can also mix the outer parts o f discs of the host 
galaxy (lOuillen et al.L 120091: iBird et al. , l20I2h and possibly ex- 
tend disc profiles ( lYounger et al.L 120071) . It would be interesting 
to study the different signatures of this mechanism compared to 
the effects of secular evolution we present here. 
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